From 0391b4c1bf23810eaed7bd5c435138a47fa499d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Sun, 10 Sep 2023 12:53:22 +0200 Subject: [PATCH] Drop support for image parameters in CameraFrame MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Stefan Fröse Co-Authored-By: Maximilian Linhoff --- docs/changes/2396.api.rst | 7 + examples/algorithms/dilate_image.py | 11 +- .../tutorials/calibrated_data_exploration.py | 13 +- examples/tutorials/ctapipe_handson.py | 42 ++--- examples/tutorials/ctapipe_overview.py | 6 +- examples/visualization/camera_display.py | 35 ++-- src/ctapipe/conftest.py | 24 --- src/ctapipe/containers.py | 41 +---- src/ctapipe/image/concentration.py | 41 ++--- src/ctapipe/image/extractor.py | 6 +- src/ctapipe/image/hillas.py | 17 +- src/ctapipe/image/image_processor.py | 47 ++---- src/ctapipe/image/tests/test_concentration.py | 15 +- src/ctapipe/image/tests/test_hillas.py | 156 ++++-------------- .../image/tests/test_image_processor.py | 52 ------ .../image/tests/test_timing_parameters.py | 19 ++- src/ctapipe/image/timing.py | 28 +--- src/ctapipe/io/hdf5eventsource.py | 9 +- src/ctapipe/io/tests/test_hdf5eventsource.py | 26 --- src/ctapipe/reco/hillas_intersection.py | 35 ++-- src/ctapipe/reco/hillas_reconstructor.py | 35 ++-- .../reco/tests/test_HillasReconstructor.py | 125 +------------- src/ctapipe/tools/process.py | 4 - src/ctapipe/tools/tests/test_process.py | 5 +- src/ctapipe/visualization/tests/test_mpl.py | 41 +---- src/ctapipe/visualization/utils.py | 8 +- 26 files changed, 195 insertions(+), 653 deletions(-) create mode 100644 docs/changes/2396.api.rst diff --git a/docs/changes/2396.api.rst b/docs/changes/2396.api.rst new file mode 100644 index 00000000000..88c156079b3 --- /dev/null +++ b/docs/changes/2396.api.rst @@ -0,0 +1,7 @@ +Drop support for computing image parameters in ``CameraFrame`` (meters on the camera plane). +Going forward, ctapipe only implements computing image parameters in telescope frame (degrees on the sky). + +* ``ctapipe-process`` and ``ImageProcessor`` no longer have options to chose between ``CameraFrame`` and ``TelescopeFrame``` +* All image parameter functions do not support inputs in length units anymore +* If using these low-level functions directly in your code, you will need to transform your camera geometries into ``TelescopeFrame`` + using ``geometry.transform_to(TelescopeFrame)``. This requires a correctly set ``frame`` attribute with a ``focal_length`` on the ``CameraGeometry``. If reading files using the ``SimTelEventSource`` or ``HDF5EventSource``, this is already the case, io plugins might need to adapt. diff --git a/examples/algorithms/dilate_image.py b/examples/algorithms/dilate_image.py index a5757be43b6..215cfde73c9 100644 --- a/examples/algorithms/dilate_image.py +++ b/examples/algorithms/dilate_image.py @@ -13,20 +13,21 @@ # %matplotlib inline from matplotlib import pyplot as plt +from ctapipe.coordinates import TelescopeFrame from ctapipe.image import dilate, tailcuts_clean, toymodel from ctapipe.instrument import SubarrayDescription from ctapipe.visualization import CameraDisplay # Load a camera from an example file subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst") -geom = subarray.tel[100].camera.geometry +geom = subarray.tel[100].camera.geometry.transform_to(TelescopeFrame()) # Create a fake camera image to display: model = toymodel.Gaussian( - x=0.2 * u.m, - y=0.0 * u.m, - width=0.05 * u.m, - length=0.15 * u.m, + x=0.2 * u.deg, + y=0.0 * u.deg, + width=0.05 * u.deg, + length=0.15 * u.deg, psi="35d", ) diff --git a/examples/tutorials/calibrated_data_exploration.py b/examples/tutorials/calibrated_data_exploration.py index 40d9d4bb202..cc3bc9ee4dd 100644 --- a/examples/tutorials/calibrated_data_exploration.py +++ b/examples/tutorials/calibrated_data_exploration.py @@ -10,6 +10,7 @@ import ctapipe from ctapipe.calib import CameraCalibrator +from ctapipe.coordinates import TelescopeFrame from ctapipe.image import hillas_parameters, tailcuts_clean from ctapipe.io import EventSource from ctapipe.utils.datasets import get_dataset_path @@ -98,7 +99,7 @@ tel_id = sorted(event.r1.tel.keys())[1] sub = source.subarray -geometry = sub.tel[tel_id].camera.geometry +geometry = sub.tel[tel_id].camera.geometry.transform_to(TelescopeFrame()) image = event.dl1.tel[tel_id].image ###################################################################### @@ -130,9 +131,8 @@ disp.overlay_moments(params, color="red", lw=3) disp.highlight_pixels(mask, color="white", alpha=0.3, linewidth=2) -plt.xlim(params.x.to_value(u.m) - 0.5, params.x.to_value(u.m) + 0.5) -plt.ylim(params.y.to_value(u.m) - 0.5, params.y.to_value(u.m) + 0.5) - +plt.xlim(params.fov_lon.to_value(u.deg) - 2, params.fov_lon.to_value(u.deg) + 2) +plt.ylim(params.fov_lat.to_value(u.deg) - 2, params.fov_lat.to_value(u.deg) + 2) ###################################################################### source.metadata @@ -166,6 +166,7 @@ cams_in_event = tels_in_event.intersection(cam_ids) first_tel_id = list(cams_in_event)[0] tel = sub.tel[first_tel_id] +tel_geom = tel.camera.geometry.transform_to(TelescopeFrame()) print("{}s in event: {}".format(tel, cams_in_event)) @@ -174,7 +175,7 @@ # image_sum = np.zeros_like( - tel.camera.geometry.pix_x.value + tel_geom.pix_x.value ) # just make an array of 0's in the same shape as the camera for tel_id in cams_in_event: @@ -187,7 +188,7 @@ plt.figure(figsize=(8, 8)) -disp = CameraDisplay(tel.camera.geometry, image=image_sum) +disp = CameraDisplay(tel_geom, image=image_sum) disp.overlay_moments(params, with_label=False) plt.title("Sum of {}x {}".format(len(cams_in_event), tel)) diff --git a/examples/tutorials/ctapipe_handson.py b/examples/tutorials/ctapipe_handson.py index 48219067be9..73b5533e517 100644 --- a/examples/tutorials/ctapipe_handson.py +++ b/examples/tutorials/ctapipe_handson.py @@ -22,6 +22,7 @@ from ctapipe import utils from ctapipe.calib import CameraCalibrator +from ctapipe.coordinates import TelescopeFrame from ctapipe.image import hillas_parameters, tailcuts_clean from ctapipe.io import EventSource, HDF5TableWriter from ctapipe.visualization import CameraDisplay @@ -129,20 +130,21 @@ def view_waveform(chan=0, pix_id=brightest_pixel): tel.optics ###################################################################### -tel.camera.geometry.pix_x +geom = tel.camera.geometry.transform_to(TelescopeFrame()) +geom.pix_x ###################################################################### -tel.camera.geometry.to_table() +geom.to_table() ###################################################################### tel.optics.mirror_area ###################################################################### -disp = CameraDisplay(tel.camera.geometry) +disp = CameraDisplay(geom) ###################################################################### -disp = CameraDisplay(tel.camera.geometry) +disp = CameraDisplay(geom) disp.image = r0tel.waveform[ 0, :, 10 ] # display channel 0, sample 0 (try others like 10) @@ -180,10 +182,10 @@ def view_waveform(chan=0, pix_id=brightest_pixel): dl1tel.peak_time ###################################################################### -CameraDisplay(tel.camera.geometry, image=dl1tel.image) +CameraDisplay(geom, image=dl1tel.image) ###################################################################### -CameraDisplay(tel.camera.geometry, image=dl1tel.peak_time) +CameraDisplay(geom, image=dl1tel.peak_time) ###################################################################### @@ -192,33 +194,33 @@ def view_waveform(chan=0, pix_id=brightest_pixel): image = dl1tel.image -mask = tailcuts_clean(tel.camera.geometry, image, picture_thresh=10, boundary_thresh=5) +mask = tailcuts_clean(geom, image, picture_thresh=10, boundary_thresh=5) mask ###################################################################### -CameraDisplay(tel.camera.geometry, image=mask) +CameraDisplay(geom, image=mask) ###################################################################### cleaned = image.copy() cleaned[~mask] = 0 ###################################################################### -disp = CameraDisplay(tel.camera.geometry, image=cleaned) +disp = CameraDisplay(geom, image=cleaned) disp.cmap = plt.cm.coolwarm disp.add_colorbar() -plt.xlim(0.5, 1.0) -plt.ylim(-1.0, 0.0) +plt.xlim(1.0, 2.0) +plt.ylim(-2.0, -1.0) ###################################################################### -params = hillas_parameters(tel.camera.geometry, cleaned) +params = hillas_parameters(geom, cleaned) print(params) ###################################################################### -disp = CameraDisplay(tel.camera.geometry, image=cleaned) +disp = CameraDisplay(geom, image=cleaned) disp.cmap = plt.cm.coolwarm disp.add_colorbar() -plt.xlim(0.5, 1.0) -plt.ylim(-1.0, 0.0) +plt.xlim(1.0, 2.0) +plt.ylim(-2.0, 1.0) disp.overlay_moments(params, color="white", lw=2) @@ -258,9 +260,10 @@ def view_waveform(chan=0, pix_id=brightest_pixel): for tel_id, tel_data in event.dl1.tel.items(): tel = source.subarray.tel[tel_id] - mask = tailcuts_clean(tel.camera.geometry, tel_data.image) + geom = tel.camera.geometry.transform_to(TelescopeFrame()) + mask = tailcuts_clean(geom, tel_data.image) if np.count_nonzero(mask) > 0: - params = hillas_parameters(tel.camera.geometry[mask], tel_data.image[mask]) + params = hillas_parameters(geom[mask], tel_data.image[mask]) ###################################################################### @@ -271,8 +274,9 @@ def view_waveform(chan=0, pix_id=brightest_pixel): for tel_id, tel_data in event.dl1.tel.items(): tel = source.subarray.tel[tel_id] - mask = tailcuts_clean(tel.camera.geometry, tel_data.image) - params = hillas_parameters(tel.camera.geometry[mask], tel_data.image[mask]) + geom = tel.camera.geometry.transform_to(TelescopeFrame()) + mask = tailcuts_clean(geom, tel_data.image) + params = hillas_parameters(geom[mask], tel_data.image[mask]) writer.write("hillas", params) diff --git a/examples/tutorials/ctapipe_overview.py b/examples/tutorials/ctapipe_overview.py index 38ed0ac3873..a52d49b3d61 100644 --- a/examples/tutorials/ctapipe_overview.py +++ b/examples/tutorials/ctapipe_overview.py @@ -24,6 +24,7 @@ from traitlets.config import Config from ctapipe.calib import CameraCalibrator +from ctapipe.coordinates import TelescopeFrame from ctapipe.image import ( ImageProcessor, camera_to_shower_coordinates, @@ -202,7 +203,7 @@ tel_id = 130 ###################################################################### -geometry = source.subarray.tel[tel_id].camera.geometry +geometry = source.subarray.tel[tel_id].camera.geometry.transform_to(TelescopeFrame()) dl1 = event.dl1.tel[tel_id] geometry, dl1 @@ -293,7 +294,7 @@ ###################################################################### long, trans = camera_to_shower_coordinates( - geometry.pix_x, geometry.pix_y, hillas.x, hillas.y, hillas.psi + geometry.pix_x, geometry.pix_y, hillas.fov_lon, hillas.fov_lat, hillas.psi ) plt.plot(long[clean], dl1.peak_time[clean], "o") @@ -539,6 +540,7 @@ ) image += new_image +geometry = geometry.transform_to(TelescopeFrame()) ###################################################################### clean = tailcuts_clean( geometry, diff --git a/examples/visualization/camera_display.py b/examples/visualization/camera_display.py index 56eece592f6..07b33e76155 100644 --- a/examples/visualization/camera_display.py +++ b/examples/visualization/camera_display.py @@ -11,7 +11,7 @@ from matplotlib.animation import FuncAnimation from matplotlib.colors import PowerNorm -from ctapipe.coordinates import CameraFrame, EngineeringCameraFrame, TelescopeFrame +from ctapipe.coordinates import EngineeringCameraFrame, TelescopeFrame from ctapipe.image import hillas_parameters, tailcuts_clean, toymodel from ctapipe.instrument import SubarrayDescription from ctapipe.io import EventSource @@ -101,7 +101,7 @@ # geom_camframe = geom -geom = geom_camframe.transform_to(EngineeringCameraFrame()) +geom = geom_camframe.transform_to(TelescopeFrame()) ###################################################################### @@ -198,9 +198,8 @@ ) disp.highlight_pixels(mask, alpha=1, linewidth=3, color="green") -ax[1].set_ylim(-0.5, 0.5) -ax[1].set_xlim(-0.5, 0.5) - +ax[1].set_ylim(-2, 2) +ax[1].set_xlim(-2, 2) ###################################################################### # Drawing a Hillas-parameter ellipse @@ -228,7 +227,7 @@ # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # # This depends on the coordinate frame of the ``CameraGeometry``. Here we -# will specify the coordinate the ``EngineerngCameraFrame``, but if you +# will specify the coordinate in the ``CameraFrame``, but if you # have enough information to do the coordinate transform, you could use # ``ICRS`` coordinates and overlay star positions. ``CameraDisplay`` will # convert the coordinate you pass in to the ``Frame`` of the display @@ -243,12 +242,8 @@ plt.figure(figsize=(6, 6)) disp = CameraDisplay(geom, image=image, cmap="gray_r") -coord = c.SkyCoord(x=0.5 * u.m, y=0.7 * u.m, frame=geom.frame) -coord_in_another_frame = c.SkyCoord(x=0.5 * u.m, y=0.7 * u.m, frame=CameraFrame()) +coord = c.SkyCoord(x=0.5 * u.m, y=0.7 * u.m, frame=geom_camframe.frame) disp.overlay_coordinate(coord, markersize=20, marker="*") -disp.overlay_coordinate( - coord_in_another_frame, markersize=20, marker="*", keep_old=True -) ###################################################################### @@ -261,11 +256,11 @@ subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst") -geom = subarray.tel[1].camera.geometry +geom = subarray.tel[1].camera.geometry.transform_to(TelescopeFrame()) -fov = 1.0 -maxwid = 0.05 -maxlen = 0.1 +fov = 2.0 +maxwid = 0.1 +maxlen = 0.5 fig, ax = plt.subplots(1, 1, figsize=(8, 6)) disp = CameraDisplay(geom, ax=ax) # we only need one display (it can be re-used) @@ -282,10 +277,10 @@ def update(frame): intens = width * length * (5e4 + 1e5 * np.random.exponential(2)) model = toymodel.Gaussian( - x=x * u.m, - y=y * u.m, - width=width * u.m, - length=length * u.m, + x=x * u.deg, + y=y * u.deg, + width=width * u.deg, + length=length * u.deg, psi=angle * u.deg, ) image, _, _ = model.generate_image( @@ -344,7 +339,7 @@ def update(frame): event = next(iter(source)) tel_id = list(event.r0.tel.keys())[0] -geom = source.subarray.tel[tel_id].camera.geometry +geom = source.subarray.tel[tel_id].camera.geometry.transform_to(TelescopeFrame()) waveform = event.r0.tel[tel_id].waveform n_chan, n_pix, n_samp = waveform.shape diff --git a/src/ctapipe/conftest.py b/src/ctapipe/conftest.py index 4d743e5d217..77f56862f0b 100644 --- a/src/ctapipe/conftest.py +++ b/src/ctapipe/conftest.py @@ -351,30 +351,6 @@ def dl1_divergent_file(dl1_tmp_path): return output -@pytest.fixture(scope="session") -def dl1_camera_frame_file(dl1_tmp_path, prod5_gamma_simtel_path): - """ - DL1 file containing both images and parameters from a gamma simulation set. - """ - from ctapipe.tools.process import ProcessorTool - - output = dl1_tmp_path / "gamma_camera_frame.dl1.h5" - - # prevent running process multiple times in case of parallel tests - with FileLock(output.with_suffix(output.suffix + ".lock")): - if output.is_file(): - return output - - argv = [ - f"--input={prod5_gamma_simtel_path}", - f"--output={output}", - "--camera-frame", - "--write-images", - ] - assert run_tool(ProcessorTool(), argv=argv, cwd=dl1_tmp_path) == 0 - return output - - @pytest.fixture(scope="session") def dl2_only_file(dl2_tmp_path, prod5_gamma_simtel_path): """ diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py index 9a5b39f4da8..3c3d0262843 100644 --- a/src/ctapipe/containers.py +++ b/src/ctapipe/containers.py @@ -33,8 +33,6 @@ "MonitoringCameraContainer", "MonitoringContainer", "MorphologyContainer", - "BaseHillasParametersContainer", - "CameraHillasParametersContainer", "CameraTimingParametersContainer", "ParticleClassificationContainer", "PedestalContainer", @@ -262,38 +260,7 @@ class TelEventIndexContainer(Container): tel_id = tel_id_field() -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. - """ - - default_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) - phi = Field(nan * u.deg, "polar coordinate of centroid", unit=u.deg) - - length = Field(nan * u.m, "standard deviation along the major-axis", unit=u.m) - length_uncertainty = Field(nan * u.m, "uncertainty of length", unit=u.m) - width = Field(nan * u.m, "standard spread along the minor-axis", unit=u.m) - width_uncertainty = Field(nan * u.m, "uncertainty of width", unit=u.m) - psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg) - - -class HillasParametersContainer(BaseHillasParametersContainer): +class HillasParametersContainer(Container): """ Hillas Parameters in a spherical system centered on the pointing position (TelescopeFrame). The cog position is given as offset in @@ -301,6 +268,11 @@ class HillasParametersContainer(BaseHillasParametersContainer): """ default_prefix = "hillas" + + intensity = Field(nan, "total intensity (size)") + skewness = Field(nan, "measure of the asymmetry") + kurtosis = Field(nan, "measure of the tailedness") + fov_lon = Field( nan * u.deg, "longitude angle in a spherical system centered on the pointing position", @@ -444,7 +416,6 @@ class ImageParametersContainer(Container): hillas = Field( default_factory=HillasParametersContainer, description="Hillas Parameters", - type=BaseHillasParametersContainer, ) timing = Field( default_factory=TimingParametersContainer, diff --git a/src/ctapipe/image/concentration.py b/src/ctapipe/image/concentration.py index b27e49e8727..fb3da0263b9 100644 --- a/src/ctapipe/image/concentration.py +++ b/src/ctapipe/image/concentration.py @@ -1,11 +1,7 @@ import astropy.units as u import numpy as np -from ..containers import ( - CameraHillasParametersContainer, - ConcentrationContainer, - HillasParametersContainer, -) +from ..containers import ConcentrationContainer from ..instrument import CameraGeometry from ..utils.quantities import all_to_value from .hillas import camera_to_shower_coordinates @@ -24,30 +20,17 @@ def concentration_parameters(geom: CameraGeometry, image, hillas_parameters): """ h = hillas_parameters - 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, - ) + 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 diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index ce772744c58..61b9ebd298f 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -1126,7 +1126,11 @@ def _apply_second_pass( # get projected distances along main image axis longitude, _ = camera_to_shower_coordinates( - camera_geometry.pix_x, camera_geometry.pix_y, hillas.x, hillas.y, hillas.psi + camera_geometry.pix_x, + camera_geometry.pix_y, + hillas.fov_lon, + hillas.fov_lon, + hillas.psi, ) # get the predicted times as a linear relation diff --git a/src/ctapipe/image/hillas.py b/src/ctapipe/image/hillas.py index 1dc66d03bdb..ea1611c91c7 100644 --- a/src/ctapipe/image/hillas.py +++ b/src/ctapipe/image/hillas.py @@ -8,7 +8,7 @@ import numpy as np from astropy.coordinates import Angle -from ..containers import CameraHillasParametersContainer, HillasParametersContainer +from ..containers import HillasParametersContainer HILLAS_ATOL = np.finfo(np.float64).eps @@ -178,21 +178,6 @@ 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( fov_lon=u.Quantity(cog_x, unit), fov_lat=u.Quantity(cog_y, unit), diff --git a/src/ctapipe/image/image_processor.py b/src/ctapipe/image/image_processor.py index bb810283a71..dae6edc45b0 100644 --- a/src/ctapipe/image/image_processor.py +++ b/src/ctapipe/image/image_processor.py @@ -1,24 +1,19 @@ """ High level image processing (ImageProcessor Component) """ - -from copy import deepcopy - import numpy as np from ctapipe.coordinates import TelescopeFrame from ..containers import ( ArrayEventContainer, - CameraHillasParametersContainer, - CameraTimingParametersContainer, ImageParametersContainer, IntensityStatisticsContainer, PeakTimeStatisticsContainer, TimingParametersContainer, ) from ..core import QualityQuery, TelescopeComponent -from ..core.traits import Bool, BoolTelescopeParameter, ComponentName, List +from ..core.traits import BoolTelescopeParameter, ComponentName, List from ..instrument import SubarrayDescription from .cleaning import ImageCleaner from .concentration import concentration_parameters @@ -41,15 +36,9 @@ kurtosis=np.float64(np.nan), ) DEFAULT_TIMING_PARAMETERS = TimingParametersContainer() -DEFAULT_TIMING_PARAMETERS_CAMFRAME = CameraTimingParametersContainer() DEFAULT_PEAKTIME_STATISTICS = PeakTimeStatisticsContainer() -DEFAULT_IMAGE_PARAMETERS_CAMFRAME = deepcopy(DEFAULT_IMAGE_PARAMETERS) -DEFAULT_IMAGE_PARAMETERS_CAMFRAME.hillas = CameraHillasParametersContainer() -DEFAULT_IMAGE_PARAMETERS_CAMFRAME.timing = CameraTimingParametersContainer() - - class ImageQualityQuery(QualityQuery): """for configuring image-wise data checks""" @@ -69,11 +58,6 @@ class ImageProcessor(TelescopeComponent): ImageCleaner, default_value="TailcutsImageCleaner" ).tag(config=True) - use_telescope_frame = Bool( - default_value=True, - help="Whether to calculate parameters in the telescope or camera frame", - ).tag(config=True) - apply_image_modifier = BoolTelescopeParameter( default_value=False, help="If true, apply ImageModifier to dl1 images" ).tag(config=True) @@ -106,16 +90,14 @@ def __init__( self.check_image = ImageQualityQuery(parent=self) - self.default_image_container = DEFAULT_IMAGE_PARAMETERS_CAMFRAME - if self.use_telescope_frame: - self.default_image_container = DEFAULT_IMAGE_PARAMETERS - 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 - } + self.default_image_container = DEFAULT_IMAGE_PARAMETERS + 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) @@ -148,11 +130,7 @@ def _parameterize_image( """ image_selected = image[signal_pixels] - 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 + geometry = self.telescope_frame_geometries[tel_id] # check if image can be parameterized: image_criteria = self.check_image(image=image_selected) @@ -189,10 +167,7 @@ def _parameterize_image( container_class=PeakTimeStatisticsContainer, ) else: - if self.use_telescope_frame: - timing = DEFAULT_TIMING_PARAMETERS - else: - timing = DEFAULT_TIMING_PARAMETERS_CAMFRAME + timing = DEFAULT_TIMING_PARAMETERS peak_time_statistics = DEFAULT_PEAKTIME_STATISTICS diff --git a/src/ctapipe/image/tests/test_concentration.py b/src/ctapipe/image/tests/test_concentration.py index b5ef419c376..d193f9ee780 100644 --- a/src/ctapipe/image/tests/test_concentration.py +++ b/src/ctapipe/image/tests/test_concentration.py @@ -1,13 +1,14 @@ import astropy.units as u import pytest +from ctapipe.coordinates import TelescopeFrame from ctapipe.image.concentration import concentration_parameters from ctapipe.image.hillas import hillas_parameters from ctapipe.image.tests.test_hillas import create_sample_image def test_concentration(prod5_lst): - geom = prod5_lst.camera.geometry + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) image, clean_mask = create_sample_image(psi="30d", geometry=geom) hillas = hillas_parameters(geom[clean_mask], image[clean_mask]) @@ -21,26 +22,26 @@ def test_concentration(prod5_lst): @pytest.mark.filterwarnings("error") def test_width_0(prod5_lst): - geom = prod5_lst.camera.geometry + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) image, clean_mask = create_sample_image(psi="30d", geometry=geom) hillas = hillas_parameters(geom[clean_mask], image[clean_mask]) - hillas.width = 0 * u.m + hillas.width = 0 * u.deg conc = concentration_parameters(geom, image, hillas) assert conc.core == 0 def test_no_pixels_near_cog(prod5_lst): - geom = prod5_lst.camera.geometry + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) image, clean_mask = create_sample_image(psi="30d", geometry=geom) hillas = hillas_parameters(geom[clean_mask], image[clean_mask]) # remove pixels close to cog from the cleaning pixels - clean_mask &= ((geom.pix_x - hillas.x) ** 2 + (geom.pix_y - hillas.y) ** 2) >= ( - 2 * geom.pixel_width**2 - ) + clean_mask &= ( + (geom.pix_x - hillas.fov_lon) ** 2 + (geom.pix_y - hillas.fov_lat) ** 2 + ) >= (2 * geom.pixel_width**2) conc = concentration_parameters(geom[clean_mask], image[clean_mask], hillas) assert conc.cog == 0 diff --git a/src/ctapipe/image/tests/test_hillas.py b/src/ctapipe/image/tests/test_hillas.py index 1b42bd736a8..75aceb9cec4 100644 --- a/src/ctapipe/image/tests/test_hillas.py +++ b/src/ctapipe/image/tests/test_hillas.py @@ -3,14 +3,11 @@ import numpy as np import pytest from astropy import units as u -from astropy.coordinates import Angle, SkyCoord +from astropy.coordinates import Angle from numpy import isclose from pytest import approx -from ctapipe.containers import ( - CameraHillasParametersContainer, - HillasParametersContainer, -) +from ctapipe.containers import HillasParametersContainer from ctapipe.coordinates import TelescopeFrame from ctapipe.image import tailcuts_clean, toymodel from ctapipe.image.hillas import HillasParameterizationError, hillas_parameters @@ -19,16 +16,16 @@ def create_sample_image( psi="-30d", - x=0.2 * u.m, - y=0.3 * u.m, - width=0.05 * u.m, - length=0.15 * u.m, + x=1 * u.deg, + y=1 * u.deg, + width=0.06 * u.deg, + length=0.3 * u.deg, intensity=1500, geometry=None, ): if geometry is None: s = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst") - geometry = s.tel[1].camera.geometry + geometry = s.tel[1].camera.geometry.transform_to(TelescopeFrame()) # make a toymodel shower model model = toymodel.Gaussian(x=x, y=y, width=width, length=length, psi=psi) @@ -64,7 +61,7 @@ def test_hillas_selected(prod5_lst): test Hillas-parameter routines on a sample image with selected values against a sample image with masked values set to zero """ - geom = prod5_lst.camera.geometry + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) image, clean_mask = create_sample_image(geometry=geom) image_zeros = image.copy() @@ -80,7 +77,7 @@ def test_hillas_selected(prod5_lst): def test_hillas_failure(prod5_lst): - geom = prod5_lst.camera.geometry + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) blank_image = np.zeros(geom.n_pixels) with pytest.raises(HillasParameterizationError): @@ -88,7 +85,7 @@ def test_hillas_failure(prod5_lst): def test_hillas_masked_array(prod5_lst): - geom = prod5_lst.camera.geometry + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) image, clean_mask = create_sample_image(psi="0d", geometry=geom) image_zeros = image.copy() @@ -102,30 +99,26 @@ def test_hillas_masked_array(prod5_lst): def test_hillas_container(prod5_lst): - geom = prod5_lst.camera.geometry + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) image, clean_mask = create_sample_image(psi="0d", geometry=geom) params = hillas_parameters(geom[clean_mask], image[clean_mask]) - assert isinstance(params, CameraHillasParametersContainer) - - geom_telescope_frame = geom.transform_to(TelescopeFrame()) - params = hillas_parameters(geom_telescope_frame[clean_mask], image[clean_mask]) assert isinstance(params, HillasParametersContainer) def test_with_toy(prod5_lst): rng = np.random.default_rng(42) - geom = prod5_lst.camera.geometry + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) - width = 0.03 * u.m - length = 0.15 * u.m - width_uncertainty = 0.00094 * u.m - length_uncertainty = 0.00465 * u.m + width = 0.06 * u.deg + length = 0.3 * u.deg + width_uncertainty = 0.0018 * u.deg + length_uncertainty = 0.0094 * u.deg intensity = 500 - xs = u.Quantity([0.5, 0.5, -0.5, -0.5], u.m) - ys = u.Quantity([0.5, -0.5, 0.5, -0.5], u.m) + xs = u.Quantity([1, 1, -1, -1], u.deg) + ys = u.Quantity([1, -1, 1, -1], u.deg) psis = Angle([-90, -45, 0, 45, 90], unit="deg") for x, y in zip(xs, ys): @@ -139,8 +132,8 @@ def test_with_toy(prod5_lst): result = hillas_parameters(geom, signal) - assert u.isclose(result.x, x, rtol=0.1) - assert u.isclose(result.y, y, rtol=0.1) + assert u.isclose(result.fov_lon, x, rtol=0.1) + assert u.isclose(result.fov_lat, y, rtol=0.1) assert u.isclose(result.width, width, rtol=0.1) assert u.isclose(result.width_uncertainty, width_uncertainty, rtol=0.4) @@ -156,14 +149,14 @@ def test_with_toy(prod5_lst): def test_skewness(prod5_lst): rng = np.random.default_rng(42) - geom = prod5_lst.camera.geometry + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) - width = 0.03 * u.m - length = 0.15 * u.m + width = 0.10 * u.deg + length = 0.3 * u.deg intensity = 2500 - xs = u.Quantity([0.5, 0.5, -0.5, -0.5], u.m) - ys = u.Quantity([0.5, -0.5, 0.5, -0.5], u.m) + xs = u.Quantity([1, 1, -1, -1], u.deg) + ys = u.Quantity([1, -1, 1, -1], u.deg) psis = Angle([-90, -45, 0, 45, 90], unit="deg") skews = [0, 0.3, 0.6] @@ -179,8 +172,8 @@ def test_skewness(prod5_lst): result = hillas_parameters(geom, signal) - assert u.isclose(result.x, x, rtol=0.1) - assert u.isclose(result.y, y, rtol=0.1) + assert u.isclose(result.fov_lon, x, rtol=0.1) + assert u.isclose(result.fov_lat, y, rtol=0.1) assert u.isclose(result.width, width, rtol=0.1) assert u.isclose(result.length, length, rtol=0.1) @@ -203,7 +196,7 @@ def test_skewness(prod5_lst): def test_straight_line_width_0(): """Test that hillas_parameters.width is 0 for a straight line of pixels""" # three pixels in a straight line - long = np.array([0, 1, 2]) * 0.01 + long = np.array([0, 1, 2]) * 1.25 trans = np.zeros(len(long)) pix_id = np.arange(len(long)) @@ -218,10 +211,11 @@ def test_straight_line_width_0(): geom = CameraGeometry( name="testcam", pix_id=pix_id, - pix_x=x * u.m, - pix_y=y * u.m, + pix_x=x * u.deg, + pix_y=y * u.deg, pix_type="hexagonal", - pix_area=np.full(len(pix_id), 1.0) * u.m**2, + pix_area=np.full(len(pix_id), 1.0) * u.deg**2, + frame=TelescopeFrame(), ) img = rng.poisson(5, size=len(long)) @@ -238,10 +232,11 @@ def test_single_pixel(): geom = CameraGeometry( name="testcam", pix_id=np.arange(9), - pix_x=x.ravel() * u.cm, - pix_y=y.ravel() * u.cm, + pix_x=x.ravel() * u.deg, + pix_y=y.ravel() * u.deg, pix_type="rectangular", - pix_area=np.full(9, 1.0) * u.cm**2, + pix_area=np.full(9, 1.0) * u.deg**2, + frame=TelescopeFrame(), ) image = np.zeros((3, 3)) @@ -253,82 +248,3 @@ def test_single_pixel(): assert hillas.length.value == 0 assert hillas.width.value == 0 assert np.isnan(hillas.psi) - - -def test_reconstruction_in_telescope_frame(prod5_lst): - """ - Compare the reconstruction in the telescope - and camera frame. - """ - np.random.seed(42) - - geom = prod5_lst.camera.geometry - telescope_frame = TelescopeFrame() - camera_frame = geom.frame - geom_nom = geom.transform_to(telescope_frame) - - width = 0.03 * u.m - length = 0.15 * u.m - intensity = 500 - - xs = u.Quantity([0.5, 0.5, -0.5, -0.5], u.m) - ys = u.Quantity([0.5, -0.5, 0.5, -0.5], u.m) - psis = Angle([-90, -45, 0, 45, 90], unit="deg") - - def distance(coord): - return np.sqrt(np.diff(coord.x) ** 2 + np.diff(coord.y) ** 2) / 2 - - def get_transformed_length(telescope_hillas, telescope_frame, camera_frame): - main_edges = u.Quantity([-telescope_hillas.length, telescope_hillas.length]) - main_lon = main_edges * np.cos(telescope_hillas.psi) + telescope_hillas.fov_lon - main_lat = main_edges * np.sin(telescope_hillas.psi) + telescope_hillas.fov_lat - cam_main_axis = SkyCoord( - fov_lon=main_lon, fov_lat=main_lat, frame=telescope_frame - ).transform_to(camera_frame) - transformed_length = distance(cam_main_axis) - return transformed_length - - def get_transformed_width(telescope_hillas, telescope_frame, camera_frame): - secondary_edges = u.Quantity([-telescope_hillas.width, telescope_hillas.width]) - secondary_lon = ( - secondary_edges * np.cos(telescope_hillas.psi) + telescope_result.fov_lon - ) - secondary_lat = ( - secondary_edges * np.sin(telescope_hillas.psi) + telescope_result.fov_lat - ) - cam_secondary_edges = SkyCoord( - fov_lon=secondary_lon, fov_lat=secondary_lat, frame=telescope_frame - ).transform_to(camera_frame) - transformed_width = distance(cam_secondary_edges) - return transformed_width - - for x, y in zip(xs, ys): - for psi in psis: - # generate a toy image - model = toymodel.Gaussian(x=x, y=y, width=width, length=length, psi=psi) - image, signal, noise = model.generate_image( - geom, intensity=intensity, nsb_level_pe=5 - ) - - telescope_result = hillas_parameters(geom_nom, signal) - camera_result = hillas_parameters(geom, signal) - assert camera_result.intensity == telescope_result.intensity - - # Compare results in both frames - transformed_cog = SkyCoord( - fov_lon=telescope_result.fov_lon, - fov_lat=telescope_result.fov_lat, - frame=telescope_frame, - ).transform_to(camera_frame) - assert u.isclose(transformed_cog.x, camera_result.x, rtol=0.01) - assert u.isclose(transformed_cog.y, camera_result.y, rtol=0.01) - - transformed_length = get_transformed_length( - telescope_result, telescope_frame, camera_frame - ) - assert u.isclose(transformed_length, camera_result.length, rtol=0.01) - - transformed_width = get_transformed_width( - telescope_result, telescope_frame, camera_frame - ) - assert u.isclose(transformed_width, camera_result.width, rtol=0.01) diff --git a/src/ctapipe/image/tests/test_image_processor.py b/src/ctapipe/image/tests/test_image_processor.py index 327b17aebf6..63c50265940 100644 --- a/src/ctapipe/image/tests/test_image_processor.py +++ b/src/ctapipe/image/tests/test_image_processor.py @@ -1,17 +1,9 @@ """ Tests for ImageProcessor functionality """ -from copy import deepcopy - -import astropy.units as u -import numpy as np from numpy import isfinite from ctapipe.calib import CameraCalibrator -from ctapipe.containers import ( - CameraHillasParametersContainer, - CameraTimingParametersContainer, -) from ctapipe.image import ImageProcessor from ctapipe.image.cleaning import MARSImageCleaner @@ -41,47 +33,3 @@ def test_image_processor(example_event, example_subarray): assert isfinite(dl1tel.parameters.peak_time_statistics.max) process_images.check_image.to_table() - - -def test_image_processor_camera_frame(example_event, example_subarray): - """ensure we get parameters in the camera frame if explicitly specified""" - event = deepcopy(example_event) - - calibrate = CameraCalibrator(subarray=example_subarray) - process_images = ImageProcessor( - subarray=example_subarray, - use_telescope_frame=False, - image_cleaner_type="MARSImageCleaner", - ) - - assert isinstance(process_images.clean, MARSImageCleaner) - - calibrate(event) - process_images(event) - - for dl1tel in event.dl1.tel.values(): - assert isfinite(dl1tel.image_mask.sum()) - assert isfinite(dl1tel.parameters.hillas.length.value) - dl1tel.parameters.hillas.length.to("meter") - assert isfinite(dl1tel.parameters.timing.slope.value) - assert isfinite(dl1tel.parameters.leakage.pixels_width_1) - assert isfinite(dl1tel.parameters.concentration.cog) - assert isfinite(dl1tel.parameters.morphology.n_pixels) - assert isfinite(dl1tel.parameters.intensity_statistics.max) - assert isfinite(dl1tel.parameters.peak_time_statistics.max) - - process_images.check_image.to_table() - - # set image to zeros to test invalid hillas parameters - # are in the correct frame - event = deepcopy(example_event) - calibrate(event) - for dl1 in event.dl1.tel.values(): - dl1.image = np.zeros_like(dl1.image) - - process_images(event) - for dl1 in event.dl1.tel.values(): - assert isinstance(dl1.parameters.hillas, CameraHillasParametersContainer) - assert isinstance(dl1.parameters.timing, CameraTimingParametersContainer) - assert np.isnan(dl1.parameters.hillas.length.value) - assert dl1.parameters.hillas.length.unit == u.m diff --git a/src/ctapipe/image/tests/test_timing_parameters.py b/src/ctapipe/image/tests/test_timing_parameters.py index a37c07b12db..9dd3bb9abb4 100644 --- a/src/ctapipe/image/tests/test_timing_parameters.py +++ b/src/ctapipe/image/tests/test_timing_parameters.py @@ -2,7 +2,8 @@ import numpy as np from numpy.testing import assert_allclose -from ctapipe.containers import CameraHillasParametersContainer +from ctapipe.containers import HillasParametersContainer +from ctapipe.coordinates import TelescopeFrame def test_psi_0(prod5_lst): @@ -16,8 +17,10 @@ def test_psi_0(prod5_lst): intercept = 1.0 deviation = 0.1 - geom = prod5_lst.camera.geometry - hillas = CameraHillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=0 * u.deg) + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) + hillas = HillasParametersContainer( + fov_lon=0 * u.deg, fov_lat=0 * u.deg, psi=0 * u.deg + ) random = np.random.default_rng(0) peak_time = intercept + grad * geom.pix_x.value @@ -45,9 +48,9 @@ def test_psi_20(prod5_lst): intercept = 1 deviation = 0.1 - geom = prod5_lst.camera.geometry + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) psi = 20 * u.deg - hillas = CameraHillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=psi) + hillas = HillasParametersContainer(fov_lon=0 * u.deg, fov_lat=0 * u.deg, psi=psi) random = np.random.default_rng(0) peak_time = intercept + grad * ( @@ -76,8 +79,10 @@ def test_ignore_negative(prod5_lst): intercept = 1.0 deviation = 0.1 - geom = prod5_lst.camera.geometry - hillas = CameraHillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=0 * u.deg) + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) + hillas = HillasParametersContainer( + fov_lon=0 * u.deg, fov_lat=0 * u.deg, psi=0 * u.deg + ) random = np.random.default_rng(0) peak_time = intercept + grad * geom.pix_x.value diff --git a/src/ctapipe/image/timing.py b/src/ctapipe/image/timing.py index ef3480a66d5..80ca2bfaecc 100644 --- a/src/ctapipe/image/timing.py +++ b/src/ctapipe/image/timing.py @@ -1,17 +1,11 @@ """ Image timing-based shower image parametrization. """ - import astropy.units as u import numpy as np from numba import njit -from ..containers import ( - CameraHillasParametersContainer, - CameraTimingParametersContainer, - HillasParametersContainer, - TimingParametersContainer, -) +from ..containers import TimingParametersContainer from ..fitting import lts_linear_regression from ..utils.quantities import all_to_value from .hillas import camera_to_shower_coordinates @@ -62,32 +56,22 @@ def timing_parameters(geom, image, peak_time, hillas_parameters, cleaning_mask=N raise ValueError("The non-masked pixels must verify signal >= 0") h = hillas_parameters - if isinstance(h, CameraHillasParametersContainer): - unit = h.x.unit - pix_x, pix_y, x, y, length, width = all_to_value( - geom.pix_x, geom.pix_y, h.x, h.y, h.length, h.width, unit=unit - ) - elif isinstance(h, HillasParametersContainer): - unit = h.fov_lon.unit - pix_x, pix_y, x, y, length, width = all_to_value( - geom.pix_x, geom.pix_y, h.fov_lon, h.fov_lat, h.length, h.width, unit=unit - ) + unit = h.fov_lon.unit + pix_x, pix_y, x, y = all_to_value( + geom.pix_x, geom.pix_y, h.fov_lon, h.fov_lat, unit=unit + ) longi, _ = camera_to_shower_coordinates( pix_x, pix_y, x, y, hillas_parameters.psi.to_value(u.rad) ) # re-fit using a robust-to-outlier algorithm - beta, error = lts_linear_regression(x=longi, y=peak_time, samples=5) + beta, _ = lts_linear_regression(x=longi, y=peak_time, samples=5) # error from lts_linear_regression is only for the used points, # recalculate for all points deviation = rmse(longi * beta[0] + beta[1], peak_time) - if unit.is_equivalent(u.m): - return CameraTimingParametersContainer( - slope=beta[0] / unit, intercept=beta[1], deviation=deviation - ) return TimingParametersContainer( slope=beta[0] / unit, intercept=beta[1], deviation=deviation ) diff --git a/src/ctapipe/io/hdf5eventsource.py b/src/ctapipe/io/hdf5eventsource.py index afbe2153d2c..16231a6fd0f 100644 --- a/src/ctapipe/io/hdf5eventsource.py +++ b/src/ctapipe/io/hdf5eventsource.py @@ -12,8 +12,6 @@ from ..containers import ( ArrayEventContainer, - CameraHillasParametersContainer, - CameraTimingParametersContainer, ConcentrationContainer, CoordinateFrameType, DispContainer, @@ -448,10 +446,9 @@ def _generator(self): timing_prefix = "timing" if self._is_hillas_in_camera_frame(): - hillas_cls = CameraHillasParametersContainer - timing_cls = CameraTimingParametersContainer - hillas_prefix = "camera_frame_hillas" - timing_prefix = "camera_frame_timing" + raise KeyError( + "Found DL1 parameters in camera frame. Please reprocess your files with a newer version of ctapipe." + ) param_readers = { table.name: self.reader.read( diff --git a/src/ctapipe/io/tests/test_hdf5eventsource.py b/src/ctapipe/io/tests/test_hdf5eventsource.py index eef8d2b0fa2..4579ca7161a 100644 --- a/src/ctapipe/io/tests/test_hdf5eventsource.py +++ b/src/ctapipe/io/tests/test_hdf5eventsource.py @@ -4,10 +4,6 @@ import numpy as np import pytest -from ctapipe.containers import ( - CameraHillasParametersContainer, - CameraTimingParametersContainer, -) from ctapipe.io import DataLevel, EventSource, HDF5EventSource @@ -231,28 +227,6 @@ def test_read_dl2(dl2_shower_geometry_file): assert impact.distance is not None -def test_dl1_camera_frame(dl1_camera_frame_file): - with HDF5EventSource(dl1_camera_frame_file) as s: - tel_id = None - for e in s: - for tel_id, dl1 in e.dl1.tel.items(): - assert isinstance( - dl1.parameters.hillas, CameraHillasParametersContainer - ) - assert isinstance( - dl1.parameters.timing, CameraTimingParametersContainer - ) - assert dl1.parameters.hillas.intensity is not None - - for tel_id, sim in e.simulation.tel.items(): - assert isinstance( - sim.true_parameters.hillas, CameraHillasParametersContainer - ) - assert sim.true_parameters.hillas.intensity is not None - - assert tel_id is not None, "did not test any events" - - def test_interpolate_pointing(dl1_mon_pointing_file): from ctapipe.io import HDF5EventSource diff --git a/src/ctapipe/reco/hillas_intersection.py b/src/ctapipe/reco/hillas_intersection.py index 603293ba883..747ee817f0a 100644 --- a/src/ctapipe/reco/hillas_intersection.py +++ b/src/ctapipe/reco/hillas_intersection.py @@ -15,13 +15,8 @@ import numpy as np from astropy.coordinates import AltAz, SkyCoord -from ..containers import ( - CameraHillasParametersContainer, - HillasParametersContainer, - ReconstructedGeometryContainer, -) +from ..containers import HillasParametersContainer, ReconstructedGeometryContainer from ..coordinates import ( - CameraFrame, MissingFrameAttributeWarning, NominalFrame, TelescopeFrame, @@ -209,24 +204,16 @@ def _predict(self, hillas_dict, array_pointing, telescopes_pointings=None): hillas_dict_mod = {} for tel_id, hillas in hillas_dict.items(): - if isinstance(hillas, CameraHillasParametersContainer): - focal_length = self.subarray.tel[tel_id].optics.equivalent_focal_length - camera_frame = CameraFrame( - telescope_pointing=telescopes_pointings[tel_id], - focal_length=focal_length, - ) - cog_coords = SkyCoord(x=hillas.x, y=hillas.y, frame=camera_frame) - cog_coords_nom = cog_coords.transform_to(nom_frame) - else: - telescope_frame = TelescopeFrame( - telescope_pointing=telescopes_pointings[tel_id] - ) - cog_coords = SkyCoord( - fov_lon=hillas.fov_lon, - fov_lat=hillas.fov_lat, - frame=telescope_frame, - ) - cog_coords_nom = cog_coords.transform_to(nom_frame) + telescope_frame = TelescopeFrame( + telescope_pointing=telescopes_pointings[tel_id] + ) + cog_coords = SkyCoord( + fov_lon=hillas.fov_lon, + fov_lat=hillas.fov_lat, + frame=telescope_frame, + ) + cog_coords_nom = cog_coords.transform_to(nom_frame) + hillas_dict_mod[tel_id] = HillasParametersContainer( fov_lon=cog_coords_nom.fov_lon, fov_lat=cog_coords_nom.fov_lat, diff --git a/src/ctapipe/reco/hillas_reconstructor.py b/src/ctapipe/reco/hillas_reconstructor.py index f3252993445..b55c9681c71 100644 --- a/src/ctapipe/reco/hillas_reconstructor.py +++ b/src/ctapipe/reco/hillas_reconstructor.py @@ -10,7 +10,7 @@ from astropy import units as u from astropy.coordinates import AltAz, Longitude, SkyCoord, cartesian_to_spherical -from ..containers import CameraHillasParametersContainer, ReconstructedGeometryContainer +from ..containers import ReconstructedGeometryContainer from ..coordinates import ( CameraFrame, MissingFrameAttributeWarning, @@ -257,8 +257,6 @@ def initialize_arrays(self, event, hillas_dict): az = np.empty(len(hillas_dict)) tel_ids = np.empty(len(hillas_dict), dtype=int) - hillas_in_camera_frame = False - for i, (tel_id, hillas) in enumerate(hillas_dict.items()): tel_ids[i] = tel_id @@ -266,15 +264,9 @@ def initialize_arrays(self, event, hillas_dict): alt[i] = pointing.altitude.to_value(u.rad) az[i] = pointing.azimuth.to_value(u.rad) - if isinstance(hillas, CameraHillasParametersContainer): - hillas_in_camera_frame = True - cog1[i] = hillas.x.to_value(u.m) - cog2[i] = hillas.y.to_value(u.m) - cam_radius[i] = self._cam_radius_m[tel_id] - else: - cog1[i] = hillas.fov_lon.to_value(u.deg) - cog2[i] = hillas.fov_lat.to_value(u.deg) - cam_radius[i] = self._cam_radius_deg[tel_id] + cog1[i] = hillas.fov_lon.to_value(u.deg) + cog2[i] = hillas.fov_lat.to_value(u.deg) + cam_radius[i] = self._cam_radius_deg[tel_id] psi[i] = hillas.psi.to_value(u.rad) weights[i] = hillas.intensity * hillas.length.value / hillas.width.value @@ -288,25 +280,18 @@ def initialize_arrays(self, event, hillas_dict): telescope_pointings = SkyCoord(alt=alt, az=az, unit=u.rad, frame=altaz) focal_length = u.Quantity(focal_length, u.m, copy=False) - camera_frame = CameraFrame( - telescope_pointing=telescope_pointings, focal_length=focal_length - ) telescope_frame = TelescopeFrame(telescope_pointing=telescope_pointings) p2_1 = cog1 + 0.1 * cam_radius * np.cos(psi) p2_2 = cog2 + 0.1 * cam_radius * np.sin(psi) - if hillas_in_camera_frame: - cog_coord = SkyCoord(x=cog1, y=cog2, unit=u.m, frame=camera_frame) - p2_coord = SkyCoord(x=p2_1, y=p2_2, unit=u.m, frame=camera_frame) - else: - p2_coord = SkyCoord( - fov_lon=p2_1, fov_lat=p2_2, unit=u.deg, frame=telescope_frame - ) - cog_coord = SkyCoord( - fov_lon=cog1, fov_lat=cog2, unit=u.deg, frame=telescope_frame - ) + p2_coord = SkyCoord( + fov_lon=p2_1, fov_lat=p2_2, unit=u.deg, frame=telescope_frame + ) + cog_coord = SkyCoord( + fov_lon=cog1, fov_lat=cog2, unit=u.deg, frame=telescope_frame + ) cog_coord = cog_coord.transform_to(altaz) p2_coord = p2_coord.transform_to(altaz) diff --git a/src/ctapipe/reco/tests/test_HillasReconstructor.py b/src/ctapipe/reco/tests/test_HillasReconstructor.py index d0f1126d47c..55d13491694 100644 --- a/src/ctapipe/reco/tests/test_HillasReconstructor.py +++ b/src/ctapipe/reco/tests/test_HillasReconstructor.py @@ -1,7 +1,6 @@ from copy import deepcopy import numpy as np -import pytest from astropy import units as u from astropy.coordinates import AltAz, SkyCoord from traitlets.config import Config @@ -10,9 +9,7 @@ from ctapipe.containers import HillasParametersContainer from ctapipe.coordinates import GroundFrame, altaz_to_righthanded_cartesian from ctapipe.image.image_processor import ImageProcessor -from ctapipe.io import SimTelEventSource from ctapipe.reco.hillas_reconstructor import HillasReconstructor -from ctapipe.utils import get_dataset_path def test_estimator_results(): @@ -125,7 +122,7 @@ def test_invalid_events(subarray_and_event_gamma_off_axis_500_gev): assert not result.is_valid -def test_reconstruction_against_simulation_telescope_frame( +def test_reconstruction_against_simulation( subarray_and_event_gamma_off_axis_500_gev, ): """Reconstruction is here done only in the TelescopeFrame, @@ -163,126 +160,6 @@ def test_reconstruction_against_simulation_telescope_frame( assert u.isclose(result.core_y, event.simulation.shower.core_y, atol=25 * u.m) -def test_reconstruction_against_simulation_camera_frame( - subarray_and_event_gamma_off_axis_500_gev, -): - """Reconstruction is here done only in the CameraFrame, - since the previous tests test already for the compatibility between - frames""" - - # 4-LST bright event already calibrated - # we'll clean it and parametrize it again in the TelescopeFrame - subarray, event = subarray_and_event_gamma_off_axis_500_gev - - # define reconstructor - calib = CameraCalibrator(subarray) - image_processor = ImageProcessor(subarray, use_telescope_frame=False) - reconstructor = HillasReconstructor(subarray) - - # Get shower geometry - calib(event) - image_processor(event) - reconstructor(event) - result = event.dl2.stereo.geometry[reconstructor.__class__.__name__] - - # get the reconstructed coordinates in the sky - reco_coord = SkyCoord(alt=result.alt, az=result.az, frame=AltAz()) - # get the simulated coordinates in the sky - true_coord = SkyCoord( - alt=event.simulation.shower.alt, az=event.simulation.shower.az, frame=AltAz() - ) - - # check that we are not more far than 0.1 degrees - assert reco_coord.separation(true_coord) < 0.1 * u.deg - - assert u.isclose(result.core_x, event.simulation.shower.core_x, atol=25 * u.m) - assert u.isclose(result.core_y, event.simulation.shower.core_y, atol=25 * u.m) - - -@pytest.mark.parametrize( - "filename", - [ - "gamma_divergent_LaPalma_baseline_20Zd_180Az_prod3_test.simtel.gz", - "gamma_LaPalma_baseline_20Zd_180Az_prod3b_test.simtel.gz", - ], -) -def test_CameraFrame_against_TelescopeFrame(filename): - input_file = get_dataset_path(filename) - # "gamma_divergent_LaPalma_baseline_20Zd_180Az_prod3_test.simtel.gz" - # ) - - source = SimTelEventSource( - input_file, max_events=10, focal_length_choice="EQUIVALENT" - ) - - # too few events survive for this test with the default quality criteria, - # use less restrictive ones - config = Config( - { - "StereoQualityQuery": { - "quality_criteria": [ - ("valid_width", "parameters.hillas.width.value > 0"), - ] - } - } - ) - - calib = CameraCalibrator(subarray=source.subarray) - reconstructor = HillasReconstructor(source.subarray, config=config) - image_processor_camera_frame = ImageProcessor( - source.subarray, use_telescope_frame=False - ) - image_processor_telescope_frame = ImageProcessor( - source.subarray, use_telescope_frame=True - ) - - reconstructed_events = 0 - - for event_telescope_frame in source: - calib(event_telescope_frame) - # make a copy of the calibrated event for the camera frame case - # later we clean and paramretrize the 2 events in the same way - # but in 2 different frames to check they return compatible results - event_camera_frame = deepcopy(event_telescope_frame) - - image_processor_telescope_frame(event_telescope_frame) - image_processor_camera_frame(event_camera_frame) - - reconstructor(event_camera_frame) - result_camera_frame = event_camera_frame.dl2.stereo.geometry[ - "HillasReconstructor" - ] - reconstructor(event_telescope_frame) - result_telescope_frame = event_telescope_frame.dl2.stereo.geometry[ - "HillasReconstructor" - ] - - assert result_camera_frame.is_valid == result_telescope_frame.is_valid - - if result_telescope_frame.is_valid: - reconstructed_events += 1 - - for field, cam in result_camera_frame.items(): - tel = getattr(result_telescope_frame, field) - - kwargs = dict(rtol=6e-3, equal_nan=True) - - if hasattr(cam, "unit"): - if cam.value == 0 or tel.value == 0: - kwargs["atol"] = 1e-6 * cam.unit - assert u.isclose( - cam, tel, **kwargs - ), f"attr {field} not matching, camera: {result_camera_frame!s} telescope: {result_telescope_frame!s}" - elif isinstance(cam, list): - assert cam == tel - else: - if cam == 0 or tel == 0: - kwargs["atol"] = 1e-6 - assert np.isclose(cam, tel, **kwargs) - - assert reconstructed_events > 0 # check that we reconstruct at least 1 event - - def test_angle(): from ctapipe.reco.hillas_reconstructor import angle diff --git a/src/ctapipe/tools/process.py b/src/ctapipe/tools/process.py index ef272346382..ac87288cec0 100644 --- a/src/ctapipe/tools/process.py +++ b/src/ctapipe/tools/process.py @@ -137,10 +137,6 @@ class ProcessorTool(Tool): "store DL1/Event/Telescope muon parameters in output", "don't store DL1/Event/Telescope muon parameters in output", ), - "camera-frame": ( - {"ImageProcessor": {"use_telescope_frame": False}}, - "Use camera frame for image parameters instead of telescope frame", - ), } classes = ( diff --git a/src/ctapipe/tools/tests/test_process.py b/src/ctapipe/tools/tests/test_process.py index 769574c71f4..84a750e34c1 100644 --- a/src/ctapipe/tools/tests/test_process.py +++ b/src/ctapipe/tools/tests/test_process.py @@ -77,7 +77,6 @@ def test_stage_1_dl1(tmp_path, dl1_image_file, dl1_parameters_file): f"--config={config}", f"--input={dl1_image_file}", f"--output={dl1b_from_dl1a_file}", - "--camera-frame", "--write-parameters", "--overwrite", ], @@ -103,8 +102,8 @@ def test_stage_1_dl1(tmp_path, dl1_image_file, dl1_parameters_file): "obs_id", "event_id", "tel_id", - "camera_frame_hillas_intensity", - "camera_frame_hillas_x", + "hillas_intensity", + "hillas_fov_lon", "concentration_cog", "leakage_pixels_width_1", ) diff --git a/src/ctapipe/visualization/tests/test_mpl.py b/src/ctapipe/visualization/tests/test_mpl.py index 70776678d11..63b40ac8273 100644 --- a/src/ctapipe/visualization/tests/test_mpl.py +++ b/src/ctapipe/visualization/tests/test_mpl.py @@ -11,10 +11,7 @@ from matplotlib import __version__ as mpl_version from ctapipe.calib.camera.calibrator import CameraCalibrator -from ctapipe.containers import ( - CameraHillasParametersContainer, - HillasParametersContainer, -) +from ctapipe.containers import HillasParametersContainer from ctapipe.coordinates.telescope_frame import TelescopeFrame from ctapipe.instrument import PixelShape, SubarrayDescription @@ -90,19 +87,6 @@ def test_camera_display_single(prod5_lst_cam, tmp_path): fig.savefig(tmp_path / "result.png") -def test_hillas_overlay_camera_frame(prod5_lst_cam, tmp_path): - from ctapipe.visualization import CameraDisplay - - fig, ax = plt.subplots() - disp = CameraDisplay(prod5_lst_cam, ax=ax) - hillas = CameraHillasParametersContainer( - x=0.1 * u.m, y=-0.1 * u.m, length=0.5 * u.m, width=0.2 * u.m, psi=90 * u.deg - ) - - disp.overlay_moments(hillas, color="w") - fig.savefig(tmp_path / "result.png") - - def test_hillas_overlay(prod5_lst_cam, tmp_path): from ctapipe.visualization import CameraDisplay @@ -208,16 +192,13 @@ def test_array_display(prod5_mst_nectarcam, reference_location): # ...with scalar color ad.set_vector_uv(np.array([1, 2, 3]) * u.m, np.array([1, 2, 3]) * u.m, c=3) - geom = prod5_mst_nectarcam.camera.geometry + geom = prod5_mst_nectarcam.camera.geometry.transform_to(TelescopeFrame()) rot_angle = 20 * u.deg - hillas = CameraHillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=rot_angle) - - # test using hillas params CameraFrame: - hillas_dict = { - 1: CameraHillasParametersContainer(length=100.0 * u.m, psi=90 * u.deg), - 2: CameraHillasParametersContainer(length=20000 * u.cm, psi="95deg"), - } + hillas = HillasParametersContainer( + fov_lon=0 * u.deg, fov_lat=0 * u.deg, psi=rot_angle + ) + # test using hillas params for divergent pointing in telescopeframe: grad = 2 intercept = 1 @@ -229,19 +210,11 @@ def test_array_display(prod5_mst_nectarcam, reference_location): cleaning_mask=np.ones(geom.n_pixels, dtype=bool), ) gradient_dict = {1: timing_rot20.slope.value, 2: timing_rot20.slope.value} + core_dict = { tel_id: dl1.parameters.core.psi for tel_id, dl1 in event.dl1.tel.items() } - ad.set_vector_hillas( - hillas_dict=hillas_dict, - core_dict=core_dict, - length=500, - time_gradient=gradient_dict, - angle_offset=0 * u.deg, - ) - ad.set_line_hillas(hillas_dict=hillas_dict, core_dict=core_dict, range=300) - # test using hillas params for divergent pointing in telescopeframe: hillas_dict = { 1: HillasParametersContainer( fov_lon=1.0 * u.deg, fov_lat=1.0 * u.deg, length=1.0 * u.deg, psi=90 * u.deg diff --git a/src/ctapipe/visualization/utils.py b/src/ctapipe/visualization/utils.py index 94988439f64..6cb28550d24 100644 --- a/src/ctapipe/visualization/utils.py +++ b/src/ctapipe/visualization/utils.py @@ -2,7 +2,7 @@ import numpy as np from astropy.coordinates import Angle -from ..containers import CameraHillasParametersContainer, HillasParametersContainer +from ..containers import HillasParametersContainer def build_hillas_overlay(hillas, unit, with_label=True, n_sigma=1): @@ -15,16 +15,12 @@ def build_hillas_overlay(hillas, unit, with_label=True, n_sigma=1): except u.UnitsError: raise ValueError("hillas must be in same frame as geometry") - # strip off any units if isinstance(hillas, HillasParametersContainer): cog_x = hillas.fov_lon.to_value(unit) cog_y = hillas.fov_lat.to_value(unit) - elif isinstance(hillas, CameraHillasParametersContainer): - cog_x = hillas.x.to_value(unit) - cog_y = hillas.y.to_value(unit) else: raise TypeError( - "hillas must be a (Camera)HillasParametersContainer" f", got: {hillas} " + "hillas must be a HillasParametersContainer" f", got: {hillas} " ) psi_rad = hillas.psi.to_value(u.rad)