fill missing rebin values with 0

This commit is contained in:
Lukas Klass 2020-04-17 09:57:13 +02:00
parent db13f0b073
commit 3271b0f7d3
6 changed files with 73 additions and 55 deletions

View File

@ -14,7 +14,7 @@ class SpectralQty:
A class to hold and work with spectral quantities A class to hold and work with spectral quantities
""" """
def __init__(self, wl: u.Quantity, qty: u.Quantity, extrapolate: bool = False): def __init__(self, wl: u.Quantity, qty: u.Quantity, fill_value: Union[bool, int, float] = 0):
""" """
Initialize a new spectral quantity Initialize a new spectral quantity
@ -25,8 +25,9 @@ class SpectralQty:
qty : Quantity qty : Quantity
The quantity values corresponding to the binned wavelengths. If the values are supplied without a unit, The quantity values corresponding to the binned wavelengths. If the values are supplied without a unit,
they are assumed to be dimensionless. they are assumed to be dimensionless.
extrapolate : bool fill_value : Union[bool, int, float]
Whether extrapolation should be allowed. If disabled, the spectrum will be truncated and a warning given. How to treat missing values. True enables extrapolation, False disables extrapolation and the spectrum will
be truncated. If a numeric value is given, the missing values will be filled with this value.
Returns Returns
------- -------
sqty : SpectralQty sqty : SpectralQty
@ -45,11 +46,11 @@ class SpectralQty:
self.qty = qty * u.dimensionless_unscaled self.qty = qty * u.dimensionless_unscaled
else: else:
error("Lengths not matching") error("Lengths not matching")
self._extrapolate = extrapolate self._fill_value = fill_value
@classmethod @classmethod
def fromFile(cls, file: str, wl_unit_default: u.Quantity = None, qty_unit_default: u.Quantity = None, def fromFile(cls, file: str, wl_unit_default: u.Quantity = None, qty_unit_default: u.Quantity = None,
extrapolate: bool = False) -> "SpectralQty": fill_value: Union[bool, int, float] = 0) -> "SpectralQty":
""" """
Initialize a new spectral quantity and read the values from a file Initialize a new spectral quantity and read the values from a file
@ -64,8 +65,9 @@ class SpectralQty:
Default unit to be used for the wavelength column if no units are provided by the file. Default unit to be used for the wavelength column if no units are provided by the file.
qty_unit_default : Quantity qty_unit_default : Quantity
Default unit to be used for the quantity column if no units are provided by the file. Default unit to be used for the quantity column if no units are provided by the file.
extrapolate : bool fill_value : Union[bool, int, float]
Whether extrapolation should be allowed. If disabled, the spectrum will be truncated and a warning given. How to treat missing values. True enables extrapolation, False disables extrapolation and the spectrum will
be truncated. If a numeric value is given, the missing values will be filled with this value.
Returns Returns
------- -------
sqty : SpectralQty sqty : SpectralQty
@ -89,7 +91,7 @@ class SpectralQty:
elif wl_unit_default is not None and qty_unit_default is not None: elif wl_unit_default is not None and qty_unit_default is not None:
data[data.colnames[0]].unit = wl_unit_default data[data.colnames[0]].unit = wl_unit_default
data[data.colnames[1]].unit = qty_unit_default data[data.colnames[1]].unit = qty_unit_default
return cls(data[data.colnames[0]].quantity, data[data.colnames[1]].quantity, extrapolate=extrapolate) return cls(data[data.colnames[0]].quantity, data[data.colnames[1]].quantity, fill_value=fill_value)
def __eq__(self, other) -> bool: def __eq__(self, other) -> bool:
""" """
@ -141,12 +143,11 @@ class SpectralQty:
# Summand is of type SpectralQty # Summand is of type SpectralQty
else: else:
if other.wl.unit.is_equivalent(self.wl.unit) and other.qty.unit.is_equivalent(self.qty.unit): if other.wl.unit.is_equivalent(self.wl.unit) and other.qty.unit.is_equivalent(self.qty.unit):
# Wavelengths are matching, just add the quantities
if len(self.wl) == len(other.wl) and (self.wl == other.wl).all(): if len(self.wl) == len(other.wl) and (self.wl == other.wl).all():
# Wavelengths are matching, just add the quantities
return SpectralQty(self.wl, self.qty + other.qty) return SpectralQty(self.wl, self.qty + other.qty)
# Wavelengths are not matching, rebinning needed
else: else:
# Rebin addend # Wavelengths are not matching, rebinning needed
other_rebinned = other.rebin(self.wl) other_rebinned = other.rebin(self.wl)
if len(self.wl) == len(other_rebinned.wl) and (self.wl == other_rebinned.wl).all(): if len(self.wl) == len(other_rebinned.wl) and (self.wl == other_rebinned.wl).all():
return SpectralQty(self.wl, self.qty + other_rebinned.qty) return SpectralQty(self.wl, self.qty + other_rebinned.qty)
@ -264,10 +265,15 @@ class SpectralQty:
if not wl.unit.is_equivalent(self.wl.unit): if not wl.unit.is_equivalent(self.wl.unit):
error("Mismatching units for rebinning: " + wl.unit + ", " + self.wl.unit) error("Mismatching units for rebinning: " + wl.unit + ", " + self.wl.unit)
if not self._extrapolate: if min(wl) < min(self.wl) or max(wl) > max(self.wl):
if min(wl) < min(self.wl) or max(wl) > max(self.wl): if isinstance(self._fill_value, bool):
logging.warning("Extrapolation disabled, bandwidth will be reduced.") if not self._fill_value:
# Remove new wavelengths where extrapolation would have been necessary logging.warning("Extrapolation disabled, bandwidth will be reduced.")
wl = [x.value for x in wl if min(self.wl) <= x <= max(self.wl)] * wl.unit # Remove new wavelengths where extrapolation would have been necessary
f = interp1d(self.wl, self.qty.value, fill_value="extrapolate") wl = [x.value for x in wl if min(self.wl) <= x <= max(self.wl)] * wl.unit
f = interp1d(self.wl, self.qty.value, fill_value="extrapolate")
else:
f = interp1d(self.wl, self.qty.value, fill_value=self._fill_value, bounds_error=False)
else:
f = interp1d(self.wl, self.qty.value)
return SpectralQty(wl, f(wl.to(self.wl.unit)) * self.qty.unit) return SpectralQty(wl, f(wl.to(self.wl.unit)) * self.qty.unit)

View File

@ -12,11 +12,12 @@ class TestAtmosphere(TestCase):
def test_calcSignal(self): def test_calcSignal(self):
self.assertEqual(self.atmosphere.calcSignal(), self.assertEqual(self.atmosphere.calcSignal(),
SpectralQty(np.arange(200, 208) << u.nm, [1.10e-15, 1.20e-15, 1.30e-15, 1.26e-15, 1.20e-15, SpectralQty(np.arange(200, 210) << u.nm,
1.12e-15, 1.02e-15, 0.9e-15] << u.W / ( np.array([1.10e-15, 1.20e-15, 1.30e-15, 1.26e-15, 1.20e-15, 1.12e-15, 1.02e-15,
u.m ** 2 * u.nm))) 0.9e-15, 0, 0]) << u.W / (u.m ** 2 * u.nm)))
def test_calcNoise(self): def test_calcNoise(self):
self.assertEqual(self.atmosphere.calcNoise(), self.assertEqual(self.atmosphere.calcNoise(),
SpectralQty(np.arange(200, 208) << u.nm, [1.1e-16, 1.2e-16, 1.3e-16, 1.4e-16, 1.5e-16, 1.6e-16, SpectralQty(np.arange(200, 208) << u.nm,
1.7e-16, 1.8e-16] << u.W / (u.m ** 2 * u.nm * u.sr))) np.array([1.1e-16, 1.2e-16, 1.3e-16, 1.4e-16, 1.5e-16, 1.6e-16, 1.7e-16,
1.8e-16]) << u.W / (u.m ** 2 * u.nm * u.sr)))

View File

@ -16,7 +16,9 @@ class TestBeamSplitter(TestCase):
def test___init__(self): def test___init__(self):
self.assertEqual(self.splitter.calcNoise(), self.assertEqual(self.splitter.calcNoise(),
SpectralQty(self.wl, [4.31413931e-96, 1.37122214e-95, 4.30844544e-95, 1.33846280e-94] << u.W / SpectralQty(self.wl, np.array([4.31413931e-96, 1.37122214e-95, 4.30844544e-95,
(u.m ** 2 * u.nm * u.sr))) 1.33846280e-94]) << u.W / (u.m ** 2 * u.nm * u.sr)))
self.assertEqual(self.splitter.calcSignal(), self.assertEqual(self.splitter.calcSignal(),
SpectralQty(self.wl, [1.20e-15, 1.30e-15, 1.40e-15, 1.35e-15] << u.W / (u.m ** 2 * u.nm))) SpectralQty(np.arange(200, 210, 1) << u.nm,
np.array([0, 1.20e-15, 1.30e-15, 1.40e-15, 1.35e-15, 0, 0, 0, 0, 0]) << u.W /
(u.m ** 2 * u.nm)))

View File

@ -15,7 +15,9 @@ class TestLens(TestCase):
def test___init__(self): def test___init__(self):
self.assertEqual(self.lens.calcNoise(), self.assertEqual(self.lens.calcNoise(),
SpectralQty(self.wl, [4.31413931e-96, 1.37122214e-95, 4.30844544e-95, 1.33846280e-94] << u.W / SpectralQty(self.wl, np.array([4.31413931e-96, 1.37122214e-95, 4.30844544e-95,
(u.m ** 2 * u.nm * u.sr))) 1.33846280e-94]) << u.W / (u.m ** 2 * u.nm * u.sr)))
self.assertEqual(self.lens.calcSignal(), self.assertEqual(self.lens.calcSignal(),
SpectralQty(self.wl, [1.20e-15, 1.30e-15, 1.40e-15, 1.35e-15] << u.W / (u.m ** 2 * u.nm))) SpectralQty(np.arange(200, 210, 1) << u.nm,
np.array([0, 1.20e-15, 1.30e-15, 1.40e-15, 1.35e-15, 0, 0, 0, 0, 0]) << u.W /
(u.m ** 2 * u.nm)))

View File

@ -15,7 +15,9 @@ class TestMirror(TestCase):
def test___init__(self): def test___init__(self):
self.assertEqual(self.mirror.calcNoise(), self.assertEqual(self.mirror.calcNoise(),
SpectralQty(self.wl, [4.31413931e-96, 1.37122214e-95, 4.30844544e-95, 1.33846280e-94] << u.W / SpectralQty(self.wl, np.array([4.31413931e-96, 1.37122214e-95, 4.30844544e-95,
(u.m ** 2 * u.nm * u.sr))) 1.33846280e-94]) << u.W / (u.m ** 2 * u.nm * u.sr)))
self.assertEqual(self.mirror.calcSignal(), self.assertEqual(self.mirror.calcSignal(),
SpectralQty(self.wl, [1.20e-15, 1.30e-15, 1.40e-15, 1.35e-15] << u.W / (u.m ** 2 * u.nm))) SpectralQty(np.arange(200, 210, 1) << u.nm, np.array([0, 1.20e-15, 1.30e-15, 1.40e-15,
1.35e-15, 0, 0, 0, 0, 0]) << u.W /
(u.m ** 2 * u.nm)))

View File

@ -28,20 +28,20 @@ class TestSpectralQty(TestCase):
np.arange(2.2, 3.0, 2e-1) << u.W / (u.m ** 2 * u.nm))) np.arange(2.2, 3.0, 2e-1) << u.W / (u.m ** 2 * u.nm)))
# SpectralQty # SpectralQty
self.assertEqual(self.sqty * SpectralQty(self.wl, np.arange(1, 5, 1) << u.m), self.assertEqual(self.sqty * SpectralQty(self.wl, np.arange(1, 5, 1) << u.m),
SpectralQty(self.wl, [1.1, 2.4, 3.9, 5.6] << u.W / (u.m * u.nm))) SpectralQty(self.wl, np.array([1.1, 2.4, 3.9, 5.6]) << u.W / (u.m * u.nm)))
self.assertEqual(SpectralQty(self.wl, np.arange(1, 5, 1) << u.m) * self.sqty, self.assertEqual(SpectralQty(self.wl, np.arange(1, 5, 1) << u.m) * self.sqty,
SpectralQty(self.wl, [1.1, 2.4, 3.9, 5.6] << u.W / (u.m * u.nm))) SpectralQty(self.wl, np.array([1.1, 2.4, 3.9, 5.6]) << u.W / (u.m * u.nm)))
# rebin without extrapolation and without reduction # rebin without extrapolation and without reduction
self.assertEqual( self.assertEqual(
self.sqty * SpectralQty(np.arange(199.5, 204.5, 1) << u.nm, np.arange(1, 6, 1) << u.m), self.sqty * SpectralQty(np.arange(199.5, 204.5, 1) << u.nm, np.arange(1, 6, 1) << u.m),
SpectralQty(self.wl, [1.65, 3.0, 4.55, 6.3] * u.W / (u.m * u.nm))) SpectralQty(self.wl, [1.65, 3.0, 4.55, 6.3] * u.W / (u.m * u.nm)))
# rebin without extrapolation and with reduction # rebin without extrapolation and with reduction
self.assertEqual( self.assertEqual(
self.sqty * SpectralQty(np.arange(200.5, 204.5, 1) << u.nm, np.arange(1, 5, 1) << u.m), self.sqty * SpectralQty(np.arange(200.5, 204.5, 1) << u.nm, np.arange(1, 5, 1) << u.m, fill_value=False),
SpectralQty(range(201, 204) << u.nm, [1.8, 3.25, 4.9] << u.W / (u.m * u.nm))) SpectralQty(np.arange(201, 204) << u.nm, np.array([1.8, 3.25, 4.9]) << u.W / (u.m * u.nm)))
# lambda # lambda
sqty_2 = lambda wl: 0.7 * u.dimensionless_unscaled self.assertEqual(self.sqty * (lambda wl: 0.7 * u.dimensionless_unscaled),
self.assertEqual(self.sqty * sqty_2, SpectralQty(self.wl, self.qty * 0.7)) SpectralQty(self.wl, self.qty * 0.7))
def test___sub__(self): def test___sub__(self):
# Quantity # Quantity
@ -58,12 +58,12 @@ class TestSpectralQty(TestCase):
SpectralQty(self.wl, [-0.4, -1.3, -2.2, -3.1] * u.W / (u.m ** 2 * u.nm))) SpectralQty(self.wl, [-0.4, -1.3, -2.2, -3.1] * u.W / (u.m ** 2 * u.nm)))
# rebin without extrapolation and with reduction # rebin without extrapolation and with reduction
self.assertEqual( self.assertEqual(
self.sqty - SpectralQty(np.arange(200.5, 204.5, 1) << u.nm, np.arange(1, 5, 1) << u.W / (u.m ** 2 * u.nm)), self.sqty - SpectralQty(np.arange(200.5, 204.5, 1) << u.nm, np.arange(1, 5, 1) << u.W / (u.m ** 2 * u.nm),
SpectralQty(range(201, 204) << u.nm, [-0.3, -1.2, -2.1] << u.W / (u.m ** 2 * u.nm))) fill_value=False),
SpectralQty(np.arange(201, 204) << u.nm, np.array([-0.3, -1.2, -2.1]) << u.W / (u.m ** 2 * u.nm)))
# lambda # lambda
sqty_2 = lambda wl: 1 * u.W / (u.m ** 2 * u.nm ** 2) * wl self.assertEqual(self.sqty - (lambda wl: 1 * u.W / (u.m ** 2 * u.nm ** 2) * wl),
self.assertEqual(self.sqty - sqty_2, SpectralQty(self.wl, np.array([-198.9, -199.8, -200.7, -201.6]) << u.W / (u.m ** 2 * u.nm)))
SpectralQty(self.wl, [-198.9, -199.8, -200.7, -201.6] << u.W / (u.m ** 2 * u.nm)))
def test___add__(self): def test___add__(self):
# Quantity # Quantity
@ -80,33 +80,38 @@ class TestSpectralQty(TestCase):
SpectralQty(self.wl, [2.6, 3.7, 4.8, 5.9] * u.W / (u.m ** 2 * u.nm))) SpectralQty(self.wl, [2.6, 3.7, 4.8, 5.9] * u.W / (u.m ** 2 * u.nm)))
# rebin without extrapolation and with reduction # rebin without extrapolation and with reduction
self.assertEqual( self.assertEqual(
self.sqty + SpectralQty(np.arange(200.5, 204.5, 1) << u.nm, np.arange(1, 5, 1) << u.W / (u.m ** 2 * u.nm)), self.sqty + SpectralQty(np.arange(200.5, 204.5, 1) << u.nm, np.arange(1, 5, 1) << u.W / (u.m ** 2 * u.nm),
SpectralQty(range(201, 204) << u.nm, [2.7, 3.8, 4.9] << u.W / (u.m ** 2 * u.nm))) fill_value=False),
SpectralQty(np.arange(201, 204) << u.nm, np.array([2.7, 3.8, 4.9]) << u.W / (u.m ** 2 * u.nm)))
# lambda # lambda
sqty_2 = lambda wl: 1 * u.W / (u.m ** 2 * u.nm ** 2) * wl self.assertEqual(self.sqty + (lambda wl: 1 * u.W / (u.m ** 2 * u.nm ** 2) * wl),
self.assertEqual(self.sqty + sqty_2, SpectralQty(self.wl, np.array([201.1, 202.2, 203.3, 204.4]) << u.W / (u.m ** 2 * u.nm)))
SpectralQty(self.wl, [201.1, 202.2, 203.3, 204.4] << u.W / (u.m ** 2 * u.nm)))
def test_rebinning(self): def test_rebinning(self):
# Test interpolation # Test interpolation
wl_new = np.arange(200.5, 210.5, 1) << u.nm wl_new = np.arange(200.5, 210.5, 1) << u.nm
sqty_res = SpectralQty(wl_new[:3], [1.15, 1.25, 1.35] << u.W / (u.m ** 2 * u.nm)) sqty_res = SpectralQty(wl_new[:3], np.array([1.15, 1.25, 1.35]) << u.W / (u.m ** 2 * u.nm))
sqty_rebin = self.sqty.rebin(wl_new) sqty_rebin = SpectralQty(self.wl, self.qty, fill_value=False).rebin(wl_new)
self.assertEqual(sqty_rebin, sqty_res) self.assertEqual(sqty_rebin, sqty_res)
# Test extrapolation # Test extrapolation
wl_new = np.arange(200.5, 210.5, 1) << u.nm wl_new = np.arange(200.5, 210.5, 1) << u.nm
sqty_res = SpectralQty(wl_new, [1.15, 1.25, 1.35, 1.45, 1.55, 1.65, 1.75, 1.85, sqty_res = SpectralQty(wl_new, np.array([1.15, 1.25, 1.35, 1.45, 1.55, 1.65, 1.75, 1.85,
1.95, 2.05] << u.W / (u.m ** 2 * u.nm)) 1.95, 2.05]) << u.W / (u.m ** 2 * u.nm))
sqty_extrapol = SpectralQty(self.wl, self.qty, extrapolate=True) sqty_rebin = SpectralQty(self.wl, self.qty, fill_value=True).rebin(wl_new)
sqty_rebin = sqty_extrapol.rebin(wl_new) self.assertEqual(sqty_rebin, sqty_res)
# Test fill value
wl_new = np.arange(200.5, 210.5, 1) << u.nm
sqty_res = SpectralQty(wl_new, np.array([1.15, 1.25, 1.35, 0, 0, 0, 0, 0, 0, 0]) << u.W / (u.m ** 2 * u.nm))
sqty_rebin = SpectralQty(self.wl, self.qty, fill_value=0).rebin(wl_new)
self.assertEqual(sqty_rebin, sqty_res) self.assertEqual(sqty_rebin, sqty_res)
# Test binning # Test binning
self.setUp() self.setUp()
wl_new = np.arange(200.5, 210, 2) << u.nm wl_new = np.arange(200.5, 210, 2) << u.nm
sqty_res = SpectralQty(wl_new[:2], [1.15, 1.35] << u.W / (u.m ** 2 * u.nm)) sqty_res = SpectralQty(wl_new[:2], np.array([1.15, 1.35]) << u.W / (u.m ** 2 * u.nm))
sqty_rebin = self.sqty.rebin(wl_new) sqty_rebin = SpectralQty(self.wl, self.qty, fill_value=False).rebin(wl_new)
self.assertEqual(sqty_rebin, sqty_res) self.assertEqual(sqty_rebin, sqty_res)
def test_fromFile(self): def test_fromFile(self):