Overburden-to-height conversion functions

[gschwefer]

Implements the overburden to height conversion functions in the atmosphere module

[maxnoe]

I wanted to make a simple plot of height vs. grammage using this function to check for the edge behavior for the cubic interpolation, but it flat out failed:

from ctapipe.io import EventSource
import numpy as np
import astropy.units as u

source = EventSource("dataset://gamma_prod5.simtel.zst")
grammage = np.linspace(0, 1200, 100) * u.g/u.cm**2
source.atmosphere_density_profile.height_from_overburden(grammage)

/home/maxnoe/Uni/CTA/ctapipe/ctapipe/atmosphere.py:307: RuntimeWarning: divide by zero encountered in log10
  self._height_interp(np.log10(overburden.to("g cm-2").value)), u.km
Traceback (most recent call last):
  File "/home/maxnoe/Uni/CTA/ctapipe/invalid_values.py", line 7, in <module>
    source.atmosphere_density_profile.height_from_overburden(grammage)
  File "/home/maxnoe/.local/conda/envs/cta-dev/lib/python3.11/site-packages/astropy/units/decorators.py", line 313, in wrapper
    return_ = wrapped_function(*func_args, **func_kwargs)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/maxnoe/Uni/CTA/ctapipe/ctapipe/atmosphere.py", line 307, in height_from_overburden
    self._height_interp(np.log10(overburden.to("g cm-2").value)), u.km
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/maxnoe/.local/conda/envs/cta-dev/lib/python3.11/site-packages/scipy/interpolate/_polyint.py", line 80, in __call__
    y = self._evaluate(x)
        ^^^^^^^^^^^^^^^^^
  File "/home/maxnoe/.local/conda/envs/cta-dev/lib/python3.11/site-packages/scipy/interpolate/_interpolate.py", line 755, in _evaluate
    below_bounds, above_bounds = self._check_bounds(x_new)
                                 ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/maxnoe/.local/conda/envs/cta-dev/lib/python3.11/site-packages/scipy/interpolate/_interpolate.py", line 784, in _check_bounds
    raise ValueError("A value ({}) in x_new is below "
ValueError: A value (-inf) in x_new is below the interpolation range's minimum value (-4.864756005337933).

[maxnoe]

The error is mostly misleading, because it reports the transformed by log-values, so the user cannot really tell what is wrong easily

[maxnoe]

Please also rebase vs. current main to get the docs build working

[maxnoe]

@gschwefer

I played around a bit locally and wanted to open a PR against your branch but accidentally pushed here directly.

Could you have a look at 51dbae7 and see if this ok with you? If not I'll revert

sorry

[gschwefer]

The commit looks alright. Should there be a test similar to the one you added for the other profiles as well, i.e. should we explicitly guard against inputing a negative height to the exponential model or is this a potentially relevant case?

[maxnoe]

@gschwefer Yes, I think the three classes should behave the same regarding invalid inputs, otherwise they are not interchangeable

[kosack]

Looks good, but we're still missing the inverse of AtmosphereProfile.line_of_site_integral where the zenith angle is taken into account (i.e. profile.distance_from_slant_depth(slant_depth, zenith_angle=60*u.deg) or whatever we want to call it), which is what you actually need for ImPACT, right?

[gschwefer]

ImPACT reconstructs slant depth, but what that is converted to as an output is up to whatever we decide it should. Currently, we would need a function that is height_from_slant_depth.

[gschwefer]

Another question about the functions currently there: Is the line_of_sight_integral function really correct? Because the distance you put into it is a distance from the observer, but the integral functions uses heights a.s.l.

[maxnoe]

It seems you ran an older version of black that formats the powers slightly differently.

Running black via pre-commit should use the same version as in the CI, e.g. copy the exact command of the CI:

$ pre-commit run --show-diff-on-failure --color=always --files $(git diff origin/main --name-only)

[gschwefer]

Yep thanks for the tip, I committed from a different ctapipe installation than normal where the pre-commit hooks are not setup properly

[maxnoe]

Hi @gschwefer Could you fix the conflicts?

[gschwefer]

I'll have a look later today

[codecov]

Codecov Report

Attention: Patch coverage is 98.55072% with 2 lines in your changes are missing coverage. Please review.

Project coverage is 92.66%. Comparing base (f8d1049) to head (0799687).
Report is 7 commits behind head on main.

Files	Patch %	Lines
src/ctapipe/atmosphere.py	97.22%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2422      +/-   ##
==========================================
+ Coverage   92.59%   92.66%   +0.06%     
==========================================
  Files         232      232              
  Lines       20034    20220     +186     
==========================================
+ Hits        18551    18736     +185     
- Misses       1483     1484       +1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[maxnoe]

docs build failure is unrelated, I am investigating

[gschwefer]

@maxnoe Looks like this tests fail is a github issue to me

[maxnoe]

Yes, seems to be something wrong / changed with the mac os base images, will check

[maxnoe]

I am a bit confused: https://github.com/actions/runner-images?tab=readme-ov-file#available-images

Says the macos-latest should be a macos 14 arm based image, but after I restarted the job, it says it's macos 12 on x86...

Maybe the rollout is not yet complete? And the failed run was running on arm / macos 14.

Edit: no, both were running on macos 12

[maxnoe]

Ok, was flaky, the latest run passed

[maxnoe]

Sonarqube is complaining about duplicated code in the tests, a large chunk is exactly equal:

    fit_reference = np.array(

        [

            [0.00 * 100000, -140.508, 1178.05, 994186, 0],

            [9.75 * 100000, -18.4377, 1265.08, 708915, 0],

            [19.00 * 100000, 0.217565, 1349.22, 636143, 0],

            [46.00 * 100000, -0.000201796, 703.745, 721128, 0],

            [106.00 * 100000, 0.000763128, 1, 1.57247e10, 0],

        ]

    )

    profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference)

Could be a module constant or a fixture

[gschwefer]

I'll have a look

[ctao-dpps-sonarqube]

Analysis Details

2 Issues

• 0 Bugs
• 0 Vulnerabilities
• 2 Code Smells

Coverage and Duplications

• 98.60% Coverage (92.70% Estimated after merge)
• 0.00% Duplicated Code (0.70% Estimated after merge)

Project ID: cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

View in SonarQube