Calculate ground frame position of lst relative to mc ref position

environment.yml

@@ -3,7 +3,7 @@ channels:
   - conda-forge
   - default
 dependencies:
-  - astropy=5
+  - astropy=5.2
   - python=3.10 # nail the python version, so conda does not try upgrading / dowgrading
   - ctapipe=0.17
   - eventio

setup.cfg

@@ -30,7 +30,7 @@ package_dir =
 python_requires = >=3.8
 zip_safe = False
 install_requires=
-    astropy~=5.0
+    astropy~=5.2
     ctapipe~=0.17.0
     protozfits~=2.0
     numpy>=1.20

src/ctapipe_io_lst/__init__.py

@@ -2,11 +2,10 @@
 """
 EventSource for LSTCam protobuf-fits.fz-files.
 """
+from ctapipe.instrument.subarray import EarthLocation
 import numpy as np
 from astropy import units as u
 from pkg_resources import resource_filename
-import os
-from os import listdir
 from ctapipe.core import Provenance
 from ctapipe.instrument import (
     ReflectorShape,
@@ -25,11 +24,13 @@
 from ctapipe.io.datalevels import DataLevel
 from ctapipe.core.traits import Bool, Float, Enum, Path
 from ctapipe.containers import (
-    CoordinateFrameType, ObservingMode, PixelStatusContainer, EventType, PointingMode, R0CameraContainer, R1CameraContainer,
+    CoordinateFrameType, PixelStatusContainer, EventType, PointingMode, R0CameraContainer, R1CameraContainer,
     SchedulingBlockContainer, ObservationBlockContainer,
 )
 from ctapipe.coordinates import CameraFrame
 
+from ctapipe_io_lst.ground_frame import ground_frame_from_earth_location
+
 from .multifiles import MultiFiles
 from .containers import LSTArrayEventContainer, LSTServiceContainer, LSTEventContainer
 from .version import __version__
@@ -44,7 +45,7 @@
     TIB_DTYPE,
 )
 from .constants import (
-    HIGH_GAIN, N_GAINS, N_PIXELS, N_SAMPLES, LST1_LOCATION,
+    HIGH_GAIN, LST_LOCATIONS, N_GAINS, N_PIXELS, N_SAMPLES, LST1_LOCATION, REFERENCE_LOCATION,
 )
 
 
@@ -254,6 +255,31 @@ class LSTEventSource(EventSource):
         )
     ).tag(config=True)
 
+    reference_position_lon = Float(
+        default_value=REFERENCE_LOCATION.lon.deg,
+        help=(
+            "Longitude of the reference location for telescope GroundFrame coordinates.",
+            " Default is the roughly area weighted average of LST-1, MAGIC-1 and MAGIC-2.",
+        )
+    ).tag(config=True)
+
+    reference_position_lat = Float(
+        default_value=REFERENCE_LOCATION.lat.deg,
+        help=(
+            "Latitude of the reference location for telescope GroundFrame coordinates.",
+            " Default is the roughly area weighted average of LST-1, MAGIC-1 and MAGIC-2.",
+        )
+    ).tag(config=True)
+
+    reference_position_height = Float(
+        default_value=REFERENCE_LOCATION.height.to_value(u.m),
+        help=(
+            "Height of the reference location for telescope GroundFrame coordinates.",
+            " Default is current MC obslevel.",
+        )
+    ).tag(config=True)
+
+
     classes = [PointingSource, EventTimeCalculator, LSTR0Corrections]
 
     def __init__(self, input_url=None, **kwargs):
@@ -290,7 +316,12 @@ def __init__(self, input_url=None, **kwargs):
         self.tel_id = self.camera_config.telescope_id
         self.run_start = Time(self.camera_config.date, format='unix')
 
-        self._subarray = self.create_subarray(self.tel_id)
+        reference_location = EarthLocation(
+            lon=self.reference_position_lon * u.deg,
+            lat=self.reference_position_lat * u.deg,
+            height=self.reference_position_height * u.m,
+        )
+        self._subarray = self.create_subarray(self.tel_id, reference_location)
         self.r0_r1_calibrator = LSTR0Corrections(
             subarray=self._subarray, parent=self
         )
@@ -368,13 +399,15 @@ def datalevels(self):
         return (DataLevel.R0, )
 
     @staticmethod
-    def create_subarray(tel_id=1):
+    def create_subarray(tel_id=1, reference_location=None):
         """
         Obtain the subarray from the EventSource
         Returns
         -------
         ctapipe.instrument.SubarrayDescription
         """
+        if reference_location is None:
+            reference_location = REFERENCE_LOCATION
 
         camera_geom = load_camera_geometry()
 
@@ -399,10 +432,11 @@ def create_subarray(tel_id=1):
 
         tel_descriptions = {tel_id: lst_tel_descr}
 
-        # put LST at 0, so that it is at the given position
-        # TODO: is that the right choice for LST + MAGIC Analysis?
-        # or should we use the same relative position and reference point as the MC?
-        tel_positions = {tel_id: [0, 0, 0] * u.m}
+        xyz = ground_frame_from_earth_location(
+            LST_LOCATIONS[tel_id],
+            reference_location,
+        ).cartesian.xyz
+        tel_positions = {tel_id: xyz}
 
         subarray = SubarrayDescription(
             name=f"LST-{tel_id} subarray",

src/ctapipe_io_lst/constants.py

@@ -24,4 +24,22 @@
 PIXEL_INDEX = np.arange(N_PIXELS)
 
 #: location of lst-1 as `~astropy.coordinates.EarthLocation`
-LST1_LOCATION = EarthLocation(lon=-17.89139 * u.deg, lat=28.76139 * u.deg, height=2184 * u.m)
+#: Taken from Abelardo's Coordinates of LST-1 & MAGIC presentation
+#: https://redmine.cta-observatory.org/attachments/65827
+LST1_LOCATION = EarthLocation(
+    lon=-17.89149701 * u.deg,
+    lat=28.76152611 * u.deg,
+    # height of central pin + distance from pin to elevation axis
+    height=2184 * u.m + 15.883 * u.m
+)
+
+#: Area averaged position of LST-1, MAGIC-1 and MAGIC-2 (using 23**2 and 17**2 m2)
+REFERENCE_LOCATION = EarthLocation(
+    lon=-17.890879 * u.deg,
+    lat=28.761579 * u.deg,
+    height=2199 * u.m,  # MC obs-level
+)
+
+LST_LOCATIONS = {
+    1: LST1_LOCATION,
+}

src/ctapipe_io_lst/ground_frame.py

@@ -0,0 +1,34 @@
+"""Conversions between ground frame and earth location
+
+This will be included in ctapipe 0.19, remove when upgrading.
+"""
+from astropy.coordinates import AltAz, ITRS, CartesianRepresentation
+from ctapipe.coordinates import GroundFrame
+
+def _altaz_to_earthlocation(altaz):
+    local_itrs = altaz.transform_to(ITRS(location=altaz.location))
+    itrs = ITRS(local_itrs.cartesian + altaz.location.get_itrs().cartesian)
+    return itrs.earth_location
+
+
+def _earthlocation_to_altaz(location, reference_location):
+    # See
+    # https://docs.astropy.org/en/stable/coordinates/common_errors.html#altaz-calculations-for-earth-based-objects
+    # for why this is necessary and we cannot just do
+    # `get_itrs().transform_to(AltAz())`
+    itrs_cart = location.get_itrs().cartesian
+    itrs_ref_cart = reference_location.get_itrs().cartesian
+    local_itrs = ITRS(itrs_cart - itrs_ref_cart, location=reference_location)
+    return local_itrs.transform_to(AltAz(location=reference_location))
+
+def ground_frame_to_earth_location(ground_frame, reference_location):
+    # in astropy, x points north, y points east, so we need a minus for y.
+    cart = CartesianRepresentation(ground_frame.x, -ground_frame.y, ground_frame.z)
+    altaz = AltAz(cart, location=reference_location)
+    return _altaz_to_earthlocation(altaz)
+
+def ground_frame_from_earth_location(location, reference_location):
+    altaz = _earthlocation_to_altaz(location, reference_location)
+    x, y, z = altaz.cartesian.xyz
+    # in astropy, x points north, y points east, so we need a minus for y.
+    return GroundFrame(x=x, y=-y, z=z)

src/ctapipe_io_lst/tests/test_lsteventsource.py

@@ -103,6 +103,14 @@ def test_subarray():
     assert source.lst_service.telescope_id == 1
     assert source.lst_service.num_modules == 265
 
+    position = source.subarray.positions[1]
+    mc_position = [-6.336, 60.405, 12.5] * u.m
+
+    # mc uses slightly different reference location and z is off
+    # so only test x/y distance
+    distance = np.linalg.norm(mc_position[:2] - position[:2])
+    assert distance < 0.5 * u.m
+
     with tempfile.NamedTemporaryFile(suffix='.h5') as f:
         subarray.to_hdf(f.name)