From eb54c6eaf9b6e706e3a58e7b4ccd797265f2f87a Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Thu, 26 Oct 2023 17:07:37 +0200 Subject: [PATCH 1/5] Check pointing coordinate frame in DispReconstructor --- ctapipe/reco/sklearn.py | 29 +++++++++++++++++++++++++---- docs/changes/2431.bugfix.rst | 2 ++ 2 files changed, 27 insertions(+), 4 deletions(-) create mode 100644 docs/changes/2431.bugfix.rst diff --git a/ctapipe/reco/sklearn.py b/ctapipe/reco/sklearn.py index 976ecb6079e..320ea5d579e 100644 --- a/ctapipe/reco/sklearn.py +++ b/ctapipe/reco/sklearn.py @@ -9,7 +9,7 @@ import astropy.units as u import joblib import numpy as np -from astropy.coordinates import AltAz +from astropy.coordinates import AltAz, SkyCoord from astropy.table import QTable, Table, vstack from astropy.utils.decorators import lazyproperty from sklearn.metrics import accuracy_score, r2_score, roc_auc_score @@ -777,13 +777,34 @@ def predict_table(self, key, table: Table) -> Dict[ReconstructionProperty, Table fov_lon = table["hillas_fov_lon"].quantity + disp * np.cos(psi) fov_lat = table["hillas_fov_lat"].quantity + disp * np.sin(psi) - # FIXME: Assume constant and parallel pointing for each run self.log.warning("Assuming constant and parallel pointing for each run") + if np.all(table["subarray_pointing_frame"] == 0): + pointing_alt = table["subarray_pointing_lat"] + pointing_az = table["subarray_pointing_lon"] + elif np.all(table["subarray_pointing_frame"] == 1): + pointing_altaz = SkyCoord( + ra=table["subarray_pointing_lon"], + dec=table["subarray_pointing_lat"], + frame="icrs", + ).transform_to(AltAz()) + pointing_alt = pointing_altaz.alt + pointing_az = pointing_altaz.az + elif np.all(table["subarray_pointing_frame"] == 2): + pointing_altaz = SkyCoord( + l=table["subarray_pointing_lon"], + b=table["subarray_pointing_lat"], + frame="galactic", + ).transform_to(AltAz()) + pointing_alt = pointing_altaz.alt + pointing_az = pointing_altaz.az + else: + raise KeyError("Unknown observation coordinate frame") + alt, az = telescope_to_horizontal( lon=fov_lon, lat=fov_lat, - pointing_alt=table["subarray_pointing_lat"], - pointing_az=table["subarray_pointing_lon"], + pointing_alt=pointing_alt, + pointing_az=pointing_az, ) altaz_result = Table( diff --git a/docs/changes/2431.bugfix.rst b/docs/changes/2431.bugfix.rst new file mode 100644 index 00000000000..92920becde9 --- /dev/null +++ b/docs/changes/2431.bugfix.rst @@ -0,0 +1,2 @@ +Check the coordinate frame in which the array pointing is given +before using it in ``DispReconstructor``. From 105c5021224fe965d2d6e7ee4f18d7a1e60c116e Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Thu, 26 Oct 2023 17:35:03 +0200 Subject: [PATCH 2/5] Use CoordinateFrameType enum --- ctapipe/reco/sklearn.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/ctapipe/reco/sklearn.py b/ctapipe/reco/sklearn.py index 320ea5d579e..70566de6020 100644 --- a/ctapipe/reco/sklearn.py +++ b/ctapipe/reco/sklearn.py @@ -22,6 +22,7 @@ from ..containers import ( ArrayEventContainer, + CoordinateFrameType, DispContainer, ParticleClassificationContainer, ReconstructedEnergyContainer, @@ -778,10 +779,10 @@ def predict_table(self, key, table: Table) -> Dict[ReconstructionProperty, Table fov_lat = table["hillas_fov_lat"].quantity + disp * np.sin(psi) self.log.warning("Assuming constant and parallel pointing for each run") - if np.all(table["subarray_pointing_frame"] == 0): + if np.all(table["subarray_pointing_frame"] is CoordinateFrameType.ALTAZ): pointing_alt = table["subarray_pointing_lat"] pointing_az = table["subarray_pointing_lon"] - elif np.all(table["subarray_pointing_frame"] == 1): + elif np.all(table["subarray_pointing_frame"] is CoordinateFrameType.ICRS): pointing_altaz = SkyCoord( ra=table["subarray_pointing_lon"], dec=table["subarray_pointing_lat"], @@ -789,7 +790,7 @@ def predict_table(self, key, table: Table) -> Dict[ReconstructionProperty, Table ).transform_to(AltAz()) pointing_alt = pointing_altaz.alt pointing_az = pointing_altaz.az - elif np.all(table["subarray_pointing_frame"] == 2): + elif np.all(table["subarray_pointing_frame"] is CoordinateFrameType.GALACTIC): pointing_altaz = SkyCoord( l=table["subarray_pointing_lon"], b=table["subarray_pointing_lat"], From b868cbdfeebe6d5b20440c55acc6fe1e61027ba7 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Thu, 26 Oct 2023 17:50:27 +0200 Subject: [PATCH 3/5] Change warning message --- ctapipe/reco/sklearn.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe/reco/sklearn.py b/ctapipe/reco/sklearn.py index 70566de6020..dc43fab87aa 100644 --- a/ctapipe/reco/sklearn.py +++ b/ctapipe/reco/sklearn.py @@ -778,7 +778,7 @@ def predict_table(self, key, table: Table) -> Dict[ReconstructionProperty, Table fov_lon = table["hillas_fov_lon"].quantity + disp * np.cos(psi) fov_lat = table["hillas_fov_lat"].quantity + disp * np.sin(psi) - self.log.warning("Assuming constant and parallel pointing for each run") + self.log.warning("FIXME: Assuming constant and parallel pointing for each run") if np.all(table["subarray_pointing_frame"] is CoordinateFrameType.ALTAZ): pointing_alt = table["subarray_pointing_lat"] pointing_az = table["subarray_pointing_lon"] From c83a9c68c3e55faf9b439dcf8700e55082d60992 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Fri, 27 Oct 2023 14:58:26 +0200 Subject: [PATCH 4/5] Just raise and error if training data not in altaz --- ctapipe/reco/sklearn.py | 32 ++++------------------- ctapipe/tools/train_disp_reconstructor.py | 8 ++++++ docs/changes/2431.bugfix.rst | 4 +-- 3 files changed, 15 insertions(+), 29 deletions(-) diff --git a/ctapipe/reco/sklearn.py b/ctapipe/reco/sklearn.py index dc43fab87aa..976ecb6079e 100644 --- a/ctapipe/reco/sklearn.py +++ b/ctapipe/reco/sklearn.py @@ -9,7 +9,7 @@ import astropy.units as u import joblib import numpy as np -from astropy.coordinates import AltAz, SkyCoord +from astropy.coordinates import AltAz from astropy.table import QTable, Table, vstack from astropy.utils.decorators import lazyproperty from sklearn.metrics import accuracy_score, r2_score, roc_auc_score @@ -22,7 +22,6 @@ from ..containers import ( ArrayEventContainer, - CoordinateFrameType, DispContainer, ParticleClassificationContainer, ReconstructedEnergyContainer, @@ -778,34 +777,13 @@ def predict_table(self, key, table: Table) -> Dict[ReconstructionProperty, Table fov_lon = table["hillas_fov_lon"].quantity + disp * np.cos(psi) fov_lat = table["hillas_fov_lat"].quantity + disp * np.sin(psi) - self.log.warning("FIXME: Assuming constant and parallel pointing for each run") - if np.all(table["subarray_pointing_frame"] is CoordinateFrameType.ALTAZ): - pointing_alt = table["subarray_pointing_lat"] - pointing_az = table["subarray_pointing_lon"] - elif np.all(table["subarray_pointing_frame"] is CoordinateFrameType.ICRS): - pointing_altaz = SkyCoord( - ra=table["subarray_pointing_lon"], - dec=table["subarray_pointing_lat"], - frame="icrs", - ).transform_to(AltAz()) - pointing_alt = pointing_altaz.alt - pointing_az = pointing_altaz.az - elif np.all(table["subarray_pointing_frame"] is CoordinateFrameType.GALACTIC): - pointing_altaz = SkyCoord( - l=table["subarray_pointing_lon"], - b=table["subarray_pointing_lat"], - frame="galactic", - ).transform_to(AltAz()) - pointing_alt = pointing_altaz.alt - pointing_az = pointing_altaz.az - else: - raise KeyError("Unknown observation coordinate frame") - + # FIXME: Assume constant and parallel pointing for each run + self.log.warning("Assuming constant and parallel pointing for each run") alt, az = telescope_to_horizontal( lon=fov_lon, lat=fov_lat, - pointing_alt=pointing_alt, - pointing_az=pointing_az, + pointing_alt=table["subarray_pointing_lat"], + pointing_az=table["subarray_pointing_lon"], ) altaz_result = Table( diff --git a/ctapipe/tools/train_disp_reconstructor.py b/ctapipe/tools/train_disp_reconstructor.py index 5953bda8038..67741e0e4e2 100644 --- a/ctapipe/tools/train_disp_reconstructor.py +++ b/ctapipe/tools/train_disp_reconstructor.py @@ -1,6 +1,7 @@ import astropy.units as u import numpy as np +from ctapipe.containers import CoordinateFrameType from ctapipe.core import Tool from ctapipe.core.traits import Bool, Int, IntTelescopeParameter, Path from ctapipe.exceptions import TooFewEvents @@ -104,6 +105,13 @@ def start(self): self.log.info("Loading events for %s", tel_type) table = self._read_table(tel_type) + if not np.all( + table["subarray_pointing_frame"] == CoordinateFrameType.ALTAZ.value + ): + raise ValueError( + "Pointing information for training data has to be provided in horizontal coordinates" + ) + self.log.info("Train models on %s events", len(table)) self.cross_validate(tel_type, table) diff --git a/docs/changes/2431.bugfix.rst b/docs/changes/2431.bugfix.rst index 92920becde9..d5da9c18960 100644 --- a/docs/changes/2431.bugfix.rst +++ b/docs/changes/2431.bugfix.rst @@ -1,2 +1,2 @@ -Check the coordinate frame in which the array pointing is given -before using it in ``DispReconstructor``. +Check that the array pointing is given in horizontal coordinates +before training a ``DispReconstructor``. From 45d7738208552cda988d3e09cfcb332840d538e5 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Fri, 27 Oct 2023 15:31:28 +0200 Subject: [PATCH 5/5] Move the check of the coordinate frame to the right place --- ctapipe/tools/train_disp_reconstructor.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/ctapipe/tools/train_disp_reconstructor.py b/ctapipe/tools/train_disp_reconstructor.py index 67741e0e4e2..7edf728d1d9 100644 --- a/ctapipe/tools/train_disp_reconstructor.py +++ b/ctapipe/tools/train_disp_reconstructor.py @@ -105,13 +105,6 @@ def start(self): self.log.info("Loading events for %s", tel_type) table = self._read_table(tel_type) - if not np.all( - table["subarray_pointing_frame"] == CoordinateFrameType.ALTAZ.value - ): - raise ValueError( - "Pointing information for training data has to be provided in horizontal coordinates" - ) - self.log.info("Train models on %s events", len(table)) self.cross_validate(tel_type, table) @@ -135,6 +128,13 @@ def _read_table(self, telescope_type): f"No events after quality query for telescope type {telescope_type}" ) + if not np.all( + table["subarray_pointing_frame"] == CoordinateFrameType.ALTAZ.value + ): + raise ValueError( + "Pointing information for training data has to be provided in horizontal coordinates" + ) + table = self.models.feature_generator(table, subarray=self.loader.subarray) table[self.models.target] = self._get_true_disp(table)