Update source to ctapipe 0.17

[maxnoe]

Tests run with the exception for the tests needing calibration files. These require new calibration files produced with ctapipe 0.17.

I intend to solve this hen/egg problem by also opening the PR in lstchain and use the two branches together to produce new calibration files with what already works.

[codecov]

Codecov Report

Base: 88.12% // Head: 88.64% // Increases project coverage by +0.52% 🎉

Coverage data is based on head (dcfd41d) compared to base (479e7e7).
Patch coverage: 94.78% of modified lines in pull request are covered.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #161      +/-   ##
==========================================
+ Coverage   88.12%   88.64%   +0.52%     
==========================================
  Files          16       18       +2     
  Lines        1507     1911     +404     
==========================================
+ Hits         1328     1694     +366     
- Misses        179      217      +38     

Impacted Files	Coverage Δ
src/ctapipe_io_lst/_dev_version/__init__.py	100.00% <ø> (ø)
src/ctapipe_io_lst/anyarray_dtypes.py	100.00% <ø> (ø)
src/ctapipe_io_lst/containers.py	100.00% <ø> (ø)
src/ctapipe_io_lst/event_time.py	94.89% <ø> (ø)
src/ctapipe_io_lst/multifiles.py	91.80% <ø> (ø)
src/ctapipe_io_lst/pointing.py	90.22% <ø> (ø)
src/ctapipe_io_lst/tests/test_multifile.py	100.00% <ø> (ø)
src/ctapipe_io_lst/tests/test_version.py	100.00% <ø> (ø)
src/ctapipe_io_lst/version.py	54.54% <ø> (ø)
src/ctapipe_io_lst/ground_frame.py	68.42% <68.42%> (ø)
... and 10 more

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[maxnoe]

@FrancaCassol @rlopezcoto @moralejo I think this is ready to be reviewed, tests are passing

[maxnoe]

Ok, the positions are now calculated to the reference position calculated from official lon/lat values using abelardo's method (area weighted average with mirror areas of 23² m² for LST and 17² m² for MAGIC-{1,2}).

The calculated positions are within 0.6 m of the monte carlo positions.

The computations require astropy 5.2 and I crosschecked them with the much more straight forward code of using pyproj, but to not require an additional dependency, I went for the astropy code here, which agrees within ~1cm with pyproj.

[maxnoe]

For future reference:

from astropy.coordinates import EarthLocation, AltAz, ITRS
from astropy.time import Time
import astropy.units as u
import matplotlib.pyplot as plt
from pyproj import CRS, Transformer
import numpy as np
from ctapipe.io import EventSource


locations = {
    "LST-1": dict(
        loc=EarthLocation(lon=-17.89149701 * u.deg, lat=28.76152611 * u.deg, height=2185 * u.m),
        rough_area=23**2 * u.m**2,
        tel_id=1,
    ),
    "MAGIC-1": dict(
        loc=EarthLocation(lon=-17.8900665372 * u.deg, lat=28.7619439944 * u.deg, height=2184.763 * u.m),
        rough_area=17**2 * u.m**2,
        tel_id=2,
    ),
    "MAGIC-2": dict(
        loc=EarthLocation(lon=-17.8905595556 * u.deg, lat=28.7613094808 * u.deg, height=2186.764 * u.m),
        rough_area=17**2 * u.m**2,
        tel_id=3,
    ),
}

s = EventSource('./test_data/mc/simtel_theta_20_az_180_gdiffuse_10evts.simtel.gz')
print(s.subarray.tel_coords)

for l in locations.values():
    l['area'] = s.subarray.tel[l['tel_id']].optics.mirror_area
    l['simtel_pos'] = s.subarray.positions[l['tel_id']].to_value(u.m)


for key in ('area', 'rough_area'):
    areas = [l[key].to_value(u.m**2) for l in locations.values()]
    ref_lon = np.average([l['loc'].lon.deg for l in locations.values()], weights=areas)
    ref_lat = np.average([l['loc'].lat.deg for l in locations.values()], weights=areas)

    print(f'Reference location ({key}): {ref_lon:.6f}°, {ref_lat:.6f}°')


    crs = CRS.from_dict(dict(proj='tmerc', ellps='WGS84', lon_0=ref_lon, lat_0=ref_lat, axis='enu', units='m'))

    proj = Transformer.from_crs(crs.geodetic_crs, crs)

    print("PROJ coordinates")
    for label, d in locations.items():
        y, x = proj.transform(d['loc'].lon.deg, d['loc'].lat.deg)
        print(f"{label:10s} {x: 7.3f} m {y: 7.3f} m")
        sx, sy, sz = d['simtel_pos']
        sy = -sy
        print(f'dx={x - sx:.2f}, {y - sy:.2f}')

    print("Astropy Coordinates")
    ref = EarthLocation(lon=ref_lon * u.deg, lat=ref_lat * u.deg, height=2199 * u.m)
    ref_cart = ref.get_itrs().cartesian

    fig, ax = plt.subplots()

    print(f'Reference location: {ref_lon:.6f}°, {ref_lat:.6f}°')

    for i, (label, d) in enumerate(locations.items()):
        cart = d['loc'].get_itrs().cartesian
        altaz = ITRS(cart - ref_cart, location=ref).transform_to(AltAz(location=ref))
        x, y, z = altaz.cartesian.xyz.to_value(u.m)

        print(f"{label:10s} {x: 7.3f} m {y: 7.3f} m")

        ax.plot(x, y, 'o', color='w', mec=f'C{i}', markeredgewidth=2, label=f'{label} (astropy)')
        sx, sy, sz = d['simtel_pos']
        sy = -sy
        print(f'dx={x - sx:.2f}, {y - sy:.2f}')
        ax.plot(sx, sy, 'o', color='k', mec=f'C{i}', markeredgewidth=2, label=f'{label} (sim-config)')


    ax.legend(ncol=3, loc='lower center', bbox_to_anchor=(0.5, 1.02))

    plt.show()

[maxnoe]

In [1]: from ctapipe.io import EventSource

In [2]: obsfile = EventSource('./test_data/real/R0/20200218/LST-1.1.Run02005.0000_first50.fits.fz')

In [3]: simfile = EventSource('./test_data/mc/simtel_theta_20_az_180_gdiffuse_10evts.simtel.gz')

In [4]: simfile.subarray.positions[1]
Out[4]: <Quantity [-6.336, 60.405, 12.5  ] m>

In [5]: obsfile.subarray.positions[1]
Out[5]: <Quantity [-5.8637701 , 60.37671721,  0.88271184] m>

[maxnoe]

The difference in height is due to this issue here: cta-observatory/lst-sim-config#43