Hillas parameters in telescope frame (#1591)

* Cache files by full url/path

* Fix part file not deleted on ctrl c

* Add container for nominal reconstruction of hillas parameters

- x,y replaced by lon and lat
- units changed to deg from m

* Add method to calculate hillas parameters in the nominal frame

- Requires pixel positions in the nominal frame instead of a geom object
- Most of the code is a duplicate of the reconstruction in the camera
frame for now

* Remove second hillas parameters function

- Decide which container to return based on pixel position unit

* Add containers for nominal frame and fix hillas-related parameters

* Use image parameters in nominal frame, bump data level version

- Save image parameters calculated in the nominal frame
- This affects the hillas parameters and the timing slope
- The coordinate transformations happen in the image processor

* Fix import

* Clean up container structure, use telescope frame for image parameters

- Stage 1 and the connected tools now calculate the hillas parameters in
the telescope frame
- Hillas and TimingParametersContainer now store values as deg. Old
containers persist as CameraHillas/TimingParametersContainer
- In order to keep a single ImageParametersContainer, base containers
have been added for hillas and timing parameters
- Failing tests adapted: x/y -> lon/lat, reconstruction algorithms use
the CameraHillasParametersContainer

* Precompute camera geometries

* Fix unit test

- Value changed slightly due to differences in the camera readout
definitions (#1588)

* Rename lon/lat to fov_lon/fov_lat

* Fix use of ImageParametersContainer

- Avoid base containers and thus lacking columns in the output file

* Address codacy complaints

* Fix default ImageParametersContainer

- Sets the Hillas/Timing parameters container in the telescope frame
as default values as opposed to the base containers

* Change datamodel version

- Replace v1.0.4 with v1.1.0 since column names have changed on top of
the units

* Use fov_lon/lat instead of x/y

* Distinguish between CameraHillasParametersContainer and HillasParametersContainer

* Adapt hillas intersection to work with hillas parameters in the
telescope and camera frame

* Remove duplicated tests of the results

- Comparing values between both reconstructions should handle all cases
- Comparing results in the telescope frame to fixed values was error prone, because the image is generated in the camera frame requiring explicit conversion every time

* Remove unused declarations and imports

* Fix datamodel version

* Fix unit test

- The test file gets generated by calling stage1 on a simtel file, so
  this has to take the current datamodel into account

* Fix comparison to nan in unit tests

* Fix reading of dl1 files

- Add v1.3 to the list of supported datamodel versions
- The renaming of the parameter containers changed the default container
  prefixes (HillasParametersContainer -> CameraHillasParametersContainer
  leads to the camera_hillas prefix), breaking backwards compatibility
- This is fixed by specifying explicit prefixes when reading <v1.3 files

* Fix tests for the dl1eventsource

- Improper comparison to np.nan lead to tests effectively being skipped
- The assertion would have failed actually, because the test file
  contains nan events (which is expected)
- The new test checks if all entries are nan instead. Due to the event
  loop in the event source, this requires collecting the values in a
  list first

* Write out correct datamodel version

* Fix unit test

* Try to fix the docs by making the base containers private

* Fix codacy complaints

* Add base containers to __all__

- This seems to be necessary for the sphinx docs

* Add option to choose camera frame instead of telescope frame

* Remove unused import

* Add missing underscore

* Fix selection of camera frame geometry

* Fix test converting to wrong unit

* Add missing variable name

* Use traits for the selection of the frame

* Fix unit test

Co-authored-by: Maximilian Nöthe <[email protected]>

ctapipe/containers.py

@@ -31,6 +31,9 @@
     "MonitoringCameraContainer",
     "MonitoringContainer",
     "MorphologyContainer",
+    "BaseHillasParametersContainer",
+    "CameraHillasParametersContainer",
+    "CameraTimingParametersContainer",
     "ParticleClassificationContainer",
     "PedestalContainer",
     "PixelStatusContainer",
@@ -46,6 +49,7 @@
     "SimulatedShowerDistribution",
     "SimulationConfigContainer",
     "TelEventIndexContainer",
+    "BaseTimingParametersContainer",
     "TimingParametersContainer",
     "TriggerContainer",
     "WaveformCalibrationContainer",
@@ -107,11 +111,25 @@ class TelEventIndexContainer(Container):
     tel_id = Field(0, "telescope identifier")
 
 
-class HillasParametersContainer(Container):
-    container_prefix = "hillas"
+class BaseHillasParametersContainer(Container):
+    """
+    Base container for hillas parameters to
+    allow the CameraHillasParametersContainer to
+    be assigned to an ImageParametersContainer as well.
+    """
 
     intensity = Field(nan, "total intensity (size)")
+    skewness = Field(nan, "measure of the asymmetry")
+    kurtosis = Field(nan, "measure of the tailedness")
+
 
+class CameraHillasParametersContainer(BaseHillasParametersContainer):
+    """
+    Hillas Parameters in the camera frame. The cog position
+    is given in meter from the camera center.
+    """
+
+    container_prefix = "camera_frame_hillas"
     x = Field(nan * u.m, "centroid x coordinate", unit=u.m)
     y = Field(nan * u.m, "centroid x coordinate", unit=u.m)
     r = Field(nan * u.m, "radial coordinate of centroid", unit=u.m)
@@ -123,8 +141,33 @@ class HillasParametersContainer(Container):
     width_uncertainty = Field(nan * u.m, "uncertainty of width", unit=u.m)
     psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg)
 
-    skewness = Field(nan, "measure of the asymmetry")
-    kurtosis = Field(nan, "measure of the tailedness")
+
+class HillasParametersContainer(BaseHillasParametersContainer):
+    """
+    Hillas Parameters in a spherical system centered on the pointing position
+    (TelescopeFrame). The cog position is given as offset in
+    longitude and latitude in degree.
+    """
+
+    container_prefix = "hillas"
+    fov_lon = Field(
+        nan * u.deg,
+        "longitude angle in a spherical system centered on the pointing position",
+        unit=u.deg,
+    )
+    fov_lat = Field(
+        nan * u.deg,
+        "latitude angle in a spherical system centered on the pointing position",
+        unit=u.deg,
+    )
+    r = Field(nan * u.deg, "radial coordinate of centroid", unit=u.deg)
+    phi = Field(nan * u.deg, "polar coordinate of centroid", unit=u.deg)
+
+    length = Field(nan * u.deg, "standard deviation along the major-axis", unit=u.deg)
+    length_uncertainty = Field(nan * u.deg, "uncertainty of length", unit=u.deg)
+    width = Field(nan * u.deg, "standard spread along the minor-axis", unit=u.deg)
+    width_uncertainty = Field(nan * u.deg, "uncertainty of width", unit=u.deg)
+    psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg)
 
 
 class LeakageContainer(Container):
@@ -167,16 +210,13 @@ class ConcentrationContainer(Container):
     pixel = Field(nan, "Percentage of photo-electrons in the brightest pixel")
 
 
-class TimingParametersContainer(Container):
+class BaseTimingParametersContainer(Container):
     """
-    Slope and Intercept of a linear regression of the arrival times
-    along the shower main axis
+    Base container for timing parameters to
+    allow the CameraTimingParametersContainer to
+    be assigned to an ImageParametersContainer as well.
     """
 
-    container_prefix = "timing"
-    slope = Field(
-        nan / u.m, "Slope of arrival times along main shower axis", unit=1 / u.m
-    )
     intercept = Field(nan, "intercept of arrival times along main shower axis")
     deviation = Field(
         nan,
@@ -185,6 +225,31 @@ class TimingParametersContainer(Container):
     )
 
 
+class CameraTimingParametersContainer(BaseTimingParametersContainer):
+    """
+    Slope and Intercept of a linear regression of the arrival times
+    along the shower main axis in the camera frame.
+    """
+
+    container_prefix = "camera_frame_timing"
+    slope = Field(
+        nan / u.m, "Slope of arrival times along main shower axis", unit=1 / u.m
+    )
+
+
+class TimingParametersContainer(BaseTimingParametersContainer):
+    """
+    Slope and Intercept of a linear regression of the arrival times
+    along the shower main axis in a
+    spherical system centered on the pointing position (TelescopeFrame)
+    """
+
+    container_prefix = "timing"
+    slope = Field(
+        nan / u.deg, "Slope of arrival times along main shower axis", unit=1 / u.deg
+    )
+
+
 class MorphologyContainer(Container):
     """ Parameters related to pixels surviving image cleaning """
 
@@ -225,8 +290,16 @@ class ImageParametersContainer(Container):
     """ Collection of image parameters """
 
     container_prefix = "params"
-    hillas = Field(HillasParametersContainer(), "Hillas Parameters")
-    timing = Field(TimingParametersContainer(), "Timing Parameters")
+    hillas = Field(
+        HillasParametersContainer(),
+        "Hillas Parameters",
+        type=BaseHillasParametersContainer,
+    )
+    timing = Field(
+        TimingParametersContainer(),
+        "Timing Parameters",
+        type=BaseTimingParametersContainer,
+    )
     leakage = Field(LeakageContainer(), "Leakage Parameters")
     concentration = Field(ConcentrationContainer(), "Concentration Parameters")
     morphology = Field(MorphologyContainer(), "Image Morphology Parameters")
@@ -260,15 +333,13 @@ class DL1CameraContainer(Container):
         dtype=np.float32,
         ndim=1,
     )
-
     image_mask = Field(
         None,
         "Boolean numpy array where True means the pixel has passed cleaning."
         " Shape: (n_pixel, )",
         dtype=np.bool_,
         ndim=1,
     )
-
     parameters = Field(None, "Image parameters", type=ImageParametersContainer)
 
 

ctapipe/image/concentration.py

@@ -1,9 +1,13 @@
 import numpy as np
 import astropy.units as u
 
-from ..instrument import CameraGeometry
-from ..containers import ConcentrationContainer
+from ..containers import (
+    ConcentrationContainer,
+    HillasParametersContainer,
+    CameraHillasParametersContainer,
+)
 from .hillas import camera_to_shower_coordinates
+from ..instrument import CameraGeometry
 from ..utils.quantities import all_to_value
 
 __all__ = ["concentration_parameters"]
@@ -20,11 +24,30 @@ def concentration_parameters(geom: CameraGeometry, image, hillas_parameters):
     """
 
     h = hillas_parameters
-    unit = h.x.unit
-
-    pix_x, pix_y, x, y, length, width, pixel_width = all_to_value(
-        geom.pix_x, geom.pix_y, h.x, h.y, h.length, h.width, geom.pixel_width, unit=unit
-    )
+    if isinstance(h, CameraHillasParametersContainer):
+        unit = h.x.unit
+        pix_x, pix_y, x, y, length, width, pixel_width = all_to_value(
+            geom.pix_x,
+            geom.pix_y,
+            h.x,
+            h.y,
+            h.length,
+            h.width,
+            geom.pixel_width,
+            unit=unit,
+        )
+    elif isinstance(h, HillasParametersContainer):
+        unit = h.fov_lon.unit
+        pix_x, pix_y, x, y, length, width, pixel_width = all_to_value(
+            geom.pix_x,
+            geom.pix_y,
+            h.fov_lon,
+            h.fov_lat,
+            h.length,
+            h.width,
+            geom.pixel_width,
+            unit=unit,
+        )
 
     delta_x = pix_x - x
     delta_y = pix_y - y

ctapipe/image/hillas.py

@@ -8,16 +8,13 @@
 import numpy as np
 from astropy.coordinates import Angle
 from astropy.units import Quantity
-from ..containers import HillasParametersContainer
+from ..containers import CameraHillasParametersContainer, HillasParametersContainer
 
 
 HILLAS_ATOL = np.finfo(np.float64).eps
 
 
-__all__ = [
-    "hillas_parameters",
-    "HillasParameterizationError",
-]
+__all__ = ["hillas_parameters", "HillasParameterizationError"]
 
 
 def camera_to_shower_coordinates(x, y, cog_x, cog_y, psi):
@@ -198,9 +195,24 @@ def hillas_parameters(geom, image):
             np.sum(((((b * A) + (a * B) + (-c * C))) ** 2.0) * image)
         ) / (2 * width)
 
+    if unit.is_equivalent(u.m):
+        return CameraHillasParametersContainer(
+            x=u.Quantity(cog_x, unit),
+            y=u.Quantity(cog_y, unit),
+            r=u.Quantity(cog_r, unit),
+            phi=Angle(cog_phi, unit=u.rad),
+            intensity=size,
+            length=u.Quantity(length, unit),
+            length_uncertainty=u.Quantity(length_uncertainty, unit),
+            width=u.Quantity(width, unit),
+            width_uncertainty=u.Quantity(width_uncertainty, unit),
+            psi=Angle(psi, unit=u.rad),
+            skewness=skewness_long,
+            kurtosis=kurtosis_long,
+        )
     return HillasParametersContainer(
-        x=u.Quantity(cog_x, unit),
-        y=u.Quantity(cog_y, unit),
+        fov_lon=u.Quantity(cog_x, unit),
+        fov_lat=u.Quantity(cog_y, unit),
         r=u.Quantity(cog_r, unit),
         phi=Angle(cog_phi, unit=u.rad),
         intensity=size,

ctapipe/image/image_processor.py

@@ -1,18 +1,18 @@
 """
 High level image processing  (ImageProcessor Component)
 """
+from ctapipe.coordinates import TelescopeFrame
 import numpy as np
 
-
 from ..containers import (
     ArrayEventContainer,
-    ImageParametersContainer,
     IntensityStatisticsContainer,
-    PeakTimeStatisticsContainer,
+    ImageParametersContainer,
     TimingParametersContainer,
+    PeakTimeStatisticsContainer,
 )
 from ..core import QualityQuery, TelescopeComponent
-from ..core.traits import List, create_class_enum_trait
+from ..core.traits import Bool, List, create_class_enum_trait
 from ..instrument import SubarrayDescription
 from . import (
     ImageCleaner,
@@ -25,6 +25,7 @@
 )
 
 
+# avoid use of base containers for unparameterized images
 DEFAULT_IMAGE_PARAMETERS = ImageParametersContainer()
 DEFAULT_TRUE_IMAGE_PARAMETERS = ImageParametersContainer()
 DEFAULT_TRUE_IMAGE_PARAMETERS.intensity_statistics = IntensityStatisticsContainer(
@@ -59,6 +60,10 @@ class ImageProcessor(TelescopeComponent):
     image_cleaner_type = create_class_enum_trait(
         base_class=ImageCleaner, default_value="TailcutsImageCleaner"
     )
+    use_telescope_frame = Bool(
+        default_value=True,
+        help="Whether to calculate parameters in the telescope or camera frame",
+    ).tag(config=True)
 
     def __init__(
         self, subarray: SubarrayDescription, config=None, parent=None, **kwargs
@@ -85,6 +90,14 @@ def __init__(
             self.image_cleaner_type, subarray=subarray, parent=self
         )
         self.check_image = ImageQualityQuery(parent=self)
+        if self.use_telescope_frame:
+            telescope_frame = TelescopeFrame()
+            self.telescope_frame_geometries = {
+                tel_id: self.subarray.tel[tel_id].camera.geometry.transform_to(
+                    telescope_frame
+                )
+                for tel_id in self.subarray.tel
+            }
 
     def __call__(self, event: ArrayEventContainer):
         self._process_telescope_event(event)
@@ -94,11 +107,11 @@ def _parameterize_image(
         tel_id,
         image,
         signal_pixels,
+        geometry,
         peak_time=None,
         default=DEFAULT_IMAGE_PARAMETERS,
     ) -> ImageParametersContainer:
         """Apply image cleaning and calculate image features
-
         Parameters
         ----------
         tel_id: int
@@ -109,15 +122,12 @@ def _parameterize_image(
             image mask
         peak_time: np.ndarray
             peak time image
-
         Returns
         -------
         ImageParametersContainer:
             cleaning mask, parameters
         """
 
-        tel = self.subarray.tel[tel_id]
-        geometry = tel.camera.geometry
         image_selected = image[signal_pixels]
 
         # check if image can be parameterized:
@@ -176,9 +186,13 @@ def _process_telescope_event(self, event):
         """
         Loop over telescopes and process the calibrated images into parameters
         """
-
         for tel_id, dl1_camera in event.dl1.tel.items():
 
+            if self.use_telescope_frame:
+                # Use the transformed geometries
+                geometry = self.telescope_frame_geometries[tel_id]
+            else:
+                geometry = self.subarray.tel[tel_id].camera.geometry
             # compute image parameters only if requested to write them
             dl1_camera.image_mask = self.clean(
                 tel_id=tel_id,
@@ -191,6 +205,7 @@ def _process_telescope_event(self, event):
                 image=dl1_camera.image,
                 signal_pixels=dl1_camera.image_mask,
                 peak_time=dl1_camera.peak_time,
+                geometry=geometry,
             )
 
             self.log.debug("params: %s", dl1_camera.parameters.as_dict(recursive=True))
@@ -205,6 +220,7 @@ def _process_telescope_event(self, event):
                     tel_id,
                     image=sim_camera.true_image,
                     signal_pixels=sim_camera.true_image > 0,
+                    geometry=geometry,
                     peak_time=None,  # true image from simulation has no peak time
                     default=DEFAULT_TRUE_IMAGE_PARAMETERS,
                 )