Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Overburden-to-height conversion functions #2422

Merged
Merged
Changes from 1 commit
Commits
Show all changes
26 commits
Select commit Hold shift + click to select a range
2cb299c
abstract function and exponential profile
gschwefer Oct 25, 2023
79a9bf3
implemented in tabulated profile
gschwefer Oct 25, 2023
9a328b8
implemented for five layer atmosphere
gschwefer Oct 26, 2023
e9dd466
circle test implemented and passing
gschwefer Oct 26, 2023
e968f87
Add changelog
gschwefer Oct 26, 2023
f4162e8
fixed five layer profile function
gschwefer Oct 27, 2023
2276663
unit test civers all five layers now
gschwefer Oct 27, 2023
45fc98c
ran black
gschwefer Oct 27, 2023
4a384fb
Improve unit handling in atmosphere module, use fill values for inter…
maxnoe Oct 27, 2023
804b74b
Edge cases and tests for exp and 5 layer profiles
gschwefer Oct 30, 2023
7f3da52
Implement function to get height from slant depth
gschwefer Nov 2, 2023
6a19c11
ran black
gschwefer Nov 2, 2023
a4fd39b
Formtted again
gschwefer Nov 2, 2023
99bf1ab
correct pre-commit hooks
gschwefer Nov 3, 2023
9c74c34
Include observer altitude in LoS integral
gschwefer Nov 3, 2023
941c8eb
Test LoS integral height and zenith dependence
gschwefer Nov 3, 2023
51f13be
Test now using global units from atmosphere module
gschwefer Nov 6, 2023
a950213
Test value of LoS integral too
gschwefer Nov 7, 2023
3c2304b
classes take Tables, work with QTables internally
gschwefer Feb 23, 2024
1067eb0
Implement slant_depth_from_height function
gschwefer Feb 23, 2024
aa1fb9c
remove superfluous dots
gschwefer Feb 29, 2024
4ced98a
Rename test functions
gschwefer Apr 4, 2024
5869a42
Implement warning for very large zenith angles
gschwefer Apr 5, 2024
8653377
Include linter changes
gschwefer Apr 5, 2024
5cd88e5
Allow for array of zenith angles in 70 deg check
gschwefer Apr 5, 2024
0799687
Implement 5 layer atmosphere as fixture
gschwefer Apr 5, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Improve unit handling in atmosphere module, use fill values for inter…
…polation
maxnoe authored and gschwefer committed Apr 5, 2024
commit 4a384fb6b08fe5337f21c303b42b03efd4cd3b89
78 changes: 42 additions & 36 deletions src/ctapipe/atmosphere.py
Original file line number Diff line number Diff line change
@@ -15,7 +15,7 @@

import numpy as np
from astropy import units as u
from astropy.table import Table
from astropy.table import QTable
from scipy.interpolate import interp1d

__all__ = [
@@ -29,6 +29,9 @@
1,
}

GRAMMAGE_UNIT = u.g / u.cm**2
DENSITY_UNIT = u.g / u.cm**3


class AtmosphereDensityProfile(abc.ABC):
"""
@@ -73,7 +76,10 @@ def height_from_overburden(self, overburden: u.Quantity) -> u.Quantity:
"""

def line_of_sight_integral(
self, distance: u.Quantity, zenith_angle=0 * u.deg, output_units=u.g / u.cm**2
self,
distance: u.Quantity,
zenith_angle=0 * u.deg,
output_units=GRAMMAGE_UNIT,
):
r"""Line-of-sight integral from the shower distance to infinity, along
the direction specified by the zenith angle. This is sometimes called
@@ -144,7 +150,7 @@ def peek(self):
return fig, axis

@classmethod
def from_table(cls, table: Table):
def from_table(cls, table: QTable):
"""return a subclass of AtmosphereDensityProfile from a serialized
table"""

@@ -191,7 +197,7 @@ class ExponentialAtmosphereDensityProfile(AtmosphereDensityProfile):
"""

scale_height: u.Quantity = 8 * u.km
scale_density: u.Quantity = 0.00125 * u.g / (u.cm**3)
scale_density: u.Quantity = 0.00125 * u.g / u.cm**3

@u.quantity_input(height=u.m)
def __call__(self, height) -> u.Quantity:
@@ -206,7 +212,7 @@ def integral(
self.scale_density * self.scale_height * np.exp(-height / self.scale_height)
)

@u.quantity_input(overburden=u.g / (u.cm * u.cm))
@u.quantity_input(overburden=u.g / u.cm**2)
def height_from_overburden(self, overburden) -> u.Quantity:
return -self.scale_height * np.log(
overburden / (self.scale_height * self.scale_density)
@@ -247,41 +253,44 @@ class TableAtmosphereDensityProfile(AtmosphereDensityProfile):
load a TableAtmosphereDensityProfile from a supported EventSource
"""

def __init__(self, table: Table):
def __init__(self, table: QTable):
"""
Parameters
----------
table: Table
table : QTable
Table with columns `height`, `density`, and `column_density`
"""

for col in ["height", "density", "column_density"]:
if col not in table.colnames:
raise ValueError(f"Missing expected column: {col} in table")

self.table = table[
(table["height"] >= 0)
& (table["density"] > 0)
& (table["column_density"] > 0)
]
valid = (
(table["height"].value >= 0)
& (table["density"].value > 0)
& (table["column_density"].value > 0)
)
self.table = QTable(table[valid], copy=False)

# interpolation is done in log-y to minimize spline wobble
log_density = np.log10(self.table["density"].to_value(DENSITY_UNIT))
log_column_density = np.log10(
self.table["column_density"].to_value(GRAMMAGE_UNIT)
)
height_km = self.table["height"].to_value(u.km)

self._density_interp = interp1d(
self.table["height"].to("km").value,
np.log10(self.table["density"].to("g cm-3").value),
interp_kwargs = dict(
kind="cubic",
bounds_error=False,
)
self._density_interp = interp1d(
height_km, log_density, fill_value=(np.nan, -np.inf), **interp_kwargs
)
self._col_density_interp = interp1d(
self.table["height"].to("km").value,
np.log10(self.table["column_density"].to("g cm-2").value),
kind="cubic",
height_km, log_column_density, fill_value=(np.nan, -np.inf), **interp_kwargs
)

self._height_interp = interp1d(
np.log10(self.table["column_density"].to("g cm-2").value),
self.table["height"].to("km").value,
kind="cubic",
log_column_density, height_km, fill_value=(np.inf, np.nan), **interp_kwargs
)

# ensure it can be read back
@@ -290,21 +299,18 @@ def __init__(self, table: Table):

@u.quantity_input(height=u.m)
def __call__(self, height) -> u.Quantity:
return u.Quantity(
10 ** self._density_interp(height.to_value(u.km)), u.g / u.cm**3
)
log_density = self._density_interp(height.to_value(u.km))
return u.Quantity(10**log_density, DENSITY_UNIT, copy=False)

@u.quantity_input(height=u.m)
def integral(self, height) -> u.Quantity:
return u.Quantity(
10 ** self._col_density_interp(height.to_value(u.km)), u.g / u.cm**2
)
log_col_density = self._col_density_interp(height.to_value(u.km))
return u.Quantity(10**log_col_density, GRAMMAGE_UNIT, copy=False)

@u.quantity_input(overburden=u.g / (u.cm * u.cm))
def height_from_overburden(self, overburden) -> u.Quantity:
return u.Quantity(
self._height_interp(np.log10(overburden.to("g cm-2").value)), u.km
)
log_overburden = np.log10(overburden.to_value(GRAMMAGE_UNIT))
return u.Quantity(self._height_interp(log_overburden), u.km, copy=False)

def __repr__(self):
return (
@@ -365,12 +371,12 @@ class FiveLayerAtmosphereDensityProfile(AtmosphereDensityProfile):
A User’s Guide", 2021, Appendix F
"""

def __init__(self, table: Table):
def __init__(self, table: QTable):
self.table = table

param_a = self.table["a"].to("g/cm2")
param_b = self.table["b"].to("g/cm2")
param_c = self.table["c"].to("km")
param_a = self.table["a"].to(GRAMMAGE_UNIT)
param_b = self.table["b"].to(GRAMMAGE_UNIT)
param_c = self.table["c"].to(u.km)

# build list of column density functions and their derivatives:
self._funcs = [
@@ -396,7 +402,7 @@ def from_array(cls, array: np.ndarray, metadata: dict = None):
if array.shape != (5, 5):
raise ValueError("expected ndarray with shape (5,5)")

table = Table(
table = QTable(
array,
names=["height", "a", "b", "c", "1/c"],
units=["cm", "g/cm2", "g/cm2", "cm", "cm-1"],
17 changes: 15 additions & 2 deletions src/ctapipe/tests/test_atmosphere.py
Original file line number Diff line number Diff line change
@@ -26,8 +26,8 @@ def get_simtel_profile_from_eventsource():
return source.atmosphere_density_profile


@pytest.fixture(scope="session", name="table_profile")
def fixture_table_profile():
@pytest.fixture(scope="session")
def table_profile():
"""a table profile for testing"""
return get_simtel_profile_from_eventsource()

@@ -179,3 +179,16 @@ def test_height_overburden_circle(table_profile):
)

assert np.allclose(circle_height_table, 47 * u.km, rtol=0.0001)


def test_out_of_range(table_profile):
with pytest.warns(RuntimeWarning, match="divide by zero"):
assert np.isposinf(table_profile.height_from_overburden(0 * u.g / u.cm**2))

assert np.isnan(table_profile.height_from_overburden(2000 * u.g / u.cm**2))

assert table_profile(150 * u.km).value == 0.0
assert np.isnan(table_profile(-0.1 * u.km))

assert table_profile.integral(150 * u.km).value == 0.0
assert np.isnan(table_profile.integral(-0.1 * u.km))