fix camera coordinates & sky spherical offset calc

Author: TjarkMiener

dl1_data_handler/reader.py

@@ -25,7 +25,7 @@
 
 from astropy import units as u
 from astropy.coordinates.earth import EarthLocation
-from astropy.coordinates import SkyCoord
+from astropy.coordinates import AltAz, SkyCoord
 from astropy.table import (
     Table,
     unique,
@@ -34,6 +34,7 @@
 )
 from astropy.time import Time
 
+from ctapipe.coordinates import CameraFrame
 from ctapipe.core import Component, QualityQuery
 from ctapipe.core.traits import (
     Bool,
@@ -51,13 +52,6 @@
 from ctapipe.io import read_table
 from dl1_data_handler.image_mapper import ImageMapper
 
-# Reference (dummy) location to insert in the SkyCoord object as the default location
-#: 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
-)
 # Reference (dummy) time to insert in the SkyCoord object as the default time
 LST_EPOCH = Time("2018-10-01T00:00:00", scale="utc")
 
@@ -227,7 +221,6 @@ def __init__(
         parent=None,
         **kwargs,
     ):
-
         super().__init__(config=config, parent=parent, **kwargs)
 
         # Register the destructor to close all open files properly
@@ -576,14 +569,11 @@ def _construct_stereo_example_identifiers(self):
         simulation_info = []
         example_identifiers = []
         for file_idx, (filename, f) in enumerate(self.files.items()):
-            if self.process_type == ProcessType.Simulation:
-                # Read simulation information for each observation
-                simulation_info.append(read_table(f, "/configuration/simulation/run"))
-                # Construct the shower simulation table
-                simshower_table = read_table(f, "/simulation/event/subarray/shower")
             # Read the trigger table.
             trigger_table = read_table(f, "/dl1/event/subarray/trigger")
             if self.process_type == ProcessType.Simulation:
+                # Construct the shower simulation table
+                simshower_table = read_table(f, "/simulation/event/subarray/shower")
                 # The shower simulation table is joined with the subarray trigger table.
                 trigger_table = join(
                     left=trigger_table,
@@ -627,7 +617,7 @@ def _construct_stereo_example_identifiers(self):
                         right=tel_pointing,
                         keys=["obs_id", "tel_id"],
                     )
-                    table_per_type = self._transform_to_spherical_offsets(
+                    table_per_type = self._transform_to_cam_coord_offsets(
                         table_per_type
                     )
                 # Apply the multiplicity cut based on the telescope type
@@ -663,6 +653,24 @@ def _multiplicity_cut_subarray(table, key_colnames):
                 events.add_column(
                     true_shower_primary_class, name="true_shower_primary_class"
                 )
+                # Read simulation information for each observation
+                simulation_info_table = read_table(f, "/configuration/simulation/run")
+                # Append the simulation information to the list of simulation information
+                simulation_info.append(simulation_info_table)
+                # Assuming min_az = max_az and min_alt = max_alt
+                fix_array_pointing = simulation_info_table.copy()
+                fix_array_pointing.keep_columns(["obs_id", "min_az", "min_alt"])
+                fix_array_pointing.rename_column("min_az", "pointing_azimuth")
+                fix_array_pointing.rename_column("min_alt", "pointing_altitude")
+                # Join the prediction table with the telescope pointing table
+                events = join(
+                    left=events,
+                    right=fix_array_pointing,
+                    keys=["obs_id"],
+                )
+                # TODO: use keep_order for astropy v7.0.0
+                events.sort(["obs_id", "event_id"])
+                events = self._transform_to_sky_spher_offsets(events)
             # Appending the events to the list of example identifiers
             example_identifiers.append(events)
 
@@ -741,9 +749,9 @@ def _transform_to_log_energy(self, table):
         table.add_column(np.log10(table["true_energy"]), name="log_true_energy")
         return table
 
-    def _transform_to_spherical_offsets(self, table) -> Table:
+    def _transform_to_cam_coord_offsets(self, table) -> Table:
         """
-        Transform Alt/Az coordinates to spherical offsets w.r.t. the telescope pointing.
+        Transform Alt/Az coordinates to camera coordinate offsets w.r.t. the telescope pointing.
 
         This method converts the Alt/Az coordinates in the provided table to spherical offsets
         w.r.t. the telescope pointing. It also calculates the angular separation between the
@@ -758,34 +766,98 @@ def _transform_to_spherical_offsets(self, table) -> Table:
         --------
         table : astropy.table.Table
             A Table with the spherical offsets and the angular separation added as new columns.
-            The telescope pointing columns are removed from the table.
+        """
+
+        tel_tables = []
+        for tel_id in self.tel_ids:
+            tel_table = table.copy()
+            tel_table.keep_columns(
+                [
+                    "obs_id",
+                    "tel_id",
+                    "true_alt",
+                    "true_az",
+                    "telescope_pointing_azimuth",
+                    "telescope_pointing_altitude",
+                ]
+            )
+            tel_table = tel_table[tel_table["tel_id"] == tel_id]
+            # Set the telescope pointing of the SkyOffsetSeparation tranform to the fix pointing
+            tel_ground_frame = self.subarray.tel_coords[
+                self.subarray.tel_ids_to_indices(tel_id)
+            ]
+            fix_tel_pointing = SkyCoord(
+                tel_table["telescope_pointing_azimuth"],
+                tel_table["telescope_pointing_altitude"],
+                location=tel_ground_frame.to_earth_location(),
+                obstime=LST_EPOCH,
+            )
+            camera_frame = CameraFrame(
+                focal_length=self.subarray.tel[tel_id].optics.equivalent_focal_length,
+                rotation=self.subarray.tel[tel_id].camera.geometry.pix_rotation,
+                telescope_pointing=fix_tel_pointing,
+            )
+            # Transform the true Alt/Az coordinates to camera coordinates
+            true_direction = SkyCoord(
+                tel_table["true_az"],
+                tel_table["true_alt"],
+                location=tel_ground_frame.to_earth_location(),
+                obstime=LST_EPOCH,
+            )
+            true_cam_position = true_direction.transform_to(camera_frame)
+            true_cam_distance = np.sqrt(
+                true_cam_position.x**2 + true_cam_position.y**2
+            )
+            tel_table.add_column(true_cam_position.x, name="cam_coord_offset_x")
+            tel_table.add_column(true_cam_position.y, name="cam_coord_offset_y")
+            tel_table.add_column(true_cam_distance, name="cam_coord_distance")
+            tel_tables.append(tel_table)
+        tel_tables = vstack(tel_tables)
+        table = join(
+            left=table,
+            right=tel_tables,
+            keys=["obs_id", "event_id", "tel_id"],
+        )
+        # TODO: use keep_order for astropy v7.0.0
+        table.sort(["obs_id", "event_id", "tel_id"])
+        return table
+
+    def _transform_to_sky_spher_offsets(self, table) -> Table:
+        """
+        Transform Alt/Az coordinates to sky spherical offsets w.r.t. the telescope pointing.
+
+        This method converts the Alt/Az coordinates in the provided table to sky spherical offsets
+        w.r.t. the telescope pointing. It also calculates the angular separation between the
+        true and telescope pointing directions.
+
+        Parameters:
+        -----------
+        table : astropy.table.Table
+            A Table containing the true Alt/Az coordinates and telescope pointing.
+
+        Returns:
+        --------
+        table : astropy.table.Table
+            A Table with the spherical offsets and the angular separation added as new columns.
         """
         # Set the telescope pointing of the SkyOffsetSeparation tranform to the fix pointing
-        fix_pointing = SkyCoord(
-            table["telescope_pointing_azimuth"],
-            table["telescope_pointing_altitude"],
-            frame="altaz",
-            location=REFERENCE_LOCATION,
+        fix_array_pointing = SkyCoord(
+            table["pointing_azimuth"],
+            table["pointing_altitude"],
+            location=self.subarray.reference_location,
             obstime=LST_EPOCH,
         )
         true_direction = SkyCoord(
             table["true_az"],
             table["true_alt"],
-            frame="altaz",
-            location=REFERENCE_LOCATION,
+            location=self.subarray.reference_location,
             obstime=LST_EPOCH,
         )
         sky_offset = fix_pointing.spherical_offsets_to(true_direction)
         angular_separation = fix_pointing.separation(true_direction)
-        table.add_column(sky_offset[0], name="spherical_offset_az")
-        table.add_column(sky_offset[1], name="spherical_offset_alt")
-        table.add_column(angular_separation, name="angular_separation")
-        table.remove_columns(
-            [
-                "telescope_pointing_azimuth",
-                "telescope_pointing_altitude",
-            ]
-        )
+        table.add_column(sky_offset[0], name="sky_offset_lon")
+        table.add_column(sky_offset[1], name="sky_offset_lat")
+        table.add_column(angular_separation, name="sky_angular_separation")
         return table
 
     def get_parameters(self, batch, dl1b_parameter_list) -> np.array:
@@ -901,7 +973,9 @@ def generate_stereo_batch(self, batch_indices) -> Table:
         batch = self._append_features(batch)
         # Add blank inputs for missing telescopes in the batch
         if self.process_type == ProcessType.Simulation:
-            batch_grouped = batch.group_by(["obs_id", "event_id", "true_shower_primary_class"])
+            batch_grouped = batch.group_by(
+                ["obs_id", "event_id", "true_shower_primary_class"]
+            )
         elif self.process_type == ProcessType.Observation:
             batch_grouped = batch.group_by(["obs_id", "event_id"])
         for group_element in batch_grouped.groups:
@@ -927,26 +1001,27 @@ def generate_stereo_batch(self, batch_indices) -> Table:
                         if "features" in group_element.colnames:
                             blank_input_row["features"] = blank_input
                         if "mono_feature_vectors" in group_element.colnames:
-                            blank_input_row["mono_feature_vectors"] = (
-                                blank_mono_feature_vectors
-                            )
+                            blank_input_row[
+                                "mono_feature_vectors"
+                            ] = blank_mono_feature_vectors
                         if "stereo_feature_vectors" in group_element.colnames:
-                            blank_input_row["stereo_feature_vectors"] = (
-                                blank_stereo_feature_vectors
-                            )
+                            blank_input_row[
+                                "stereo_feature_vectors"
+                            ] = blank_stereo_feature_vectors
                         batch.add_row(blank_input_row)
         # Sort the batch with the new rows of blank inputs
         batch.sort(["obs_id", "event_id", "tel_type_id", "tel_id"])
         return batch
 
     def __destructor(self):
         """Destructor to ensure all opened HDF5 files are properly closed."""
-        if hasattr(self, "files"):  # Ensure self.files exists before attempting to close
+        if hasattr(
+            self, "files"
+        ):  # Ensure self.files exists before attempting to close
             for file_name in list(self.files.keys()):
                 if self.files[file_name].isopen:  # Check if file is still open
                     self.files[file_name].close()
 
-
     @abstractmethod
     def _append_features(self, batch) -> Table:
         pass
@@ -1078,7 +1153,6 @@ def __init__(
         parent=None,
         **kwargs,
     ):
-
         super().__init__(
             input_url_signal=input_url_signal,
             input_url_background=input_url_background,
@@ -1337,7 +1411,6 @@ def __init__(
         parent=None,
         **kwargs,
     ):
-
         super().__init__(
             input_url_signal=input_url_signal,
             input_url_background=input_url_background,