From 473240f6ee4a32460f8f78d9693fe4920a7a2767 Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Wed, 9 Jan 2019 18:50:54 +0100 Subject: [PATCH 1/7] deleted summary.html - accidentally committed file (#919) --- summary.html | 1019 -------------------------------------------------- 1 file changed, 1019 deletions(-) delete mode 100644 summary.html diff --git a/summary.html b/summary.html deleted file mode 100644 index 4777d102f4d..00000000000 --- a/summary.html +++ /dev/null @@ -1,1019 +0,0 @@ - -Run Summary - -

Summary for runs 127755 to 135869

Dates: 2017-01-16 to 2018-01-05

Hours per Target

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Stereo 4-telStereo 3-telMono CT5Hybrid 5-telHybrid 4-telHybrid 3-telTOTAL/TARG
1RXS J002200.9+0006596.21.47.6
1RXS J012338.2-2311006.06.0
1RXS J023832.6-3116580.34.78.81.515.3
1RXS J054357.3-5532061.93.80.96.6
3C 2790.911.312.2
3C2791.91.9
3C279ToO0.50.5
3FGL J0244.8-581819.114.934.1
3FGL J0536.4-33474.212.616.8
3FGL J0812.0+02373.12.05.1
3FGL J1936.9-47190.98.89.7
3FGL J2131.5-09156.121.127.1
ASASSN-17mt8.411.319.7
CAL 601.919.921.7
CAL600.92.63.6
Crab Nebula0.45.45.8
Crab Nebula off0.50.40.4
Crab Nebula off1.02.12.1
Crab Nebula off1.250.91.92.8
Crab Nebula off1.51.91.9
Crab Nebula off2.00.90.9
Crab nebula1.94.46.2
Eta Carinae0.227.928.1
Fornax A0.50.5
Fornax A west lobe21.816.48.846.9
G2975955.35.3
G2980484.24.2
GRB 170424A0.80.8
GRB1704020.40.4
GRB1705311.91.9
GRB1708261.91.9
GRB171112A3.13.1
GRB1712231.91.9
GW1708171.43.14.5
H 1914-1941.10.913.61.417.0
HESS J0632+0570.80.45.06.2
HESS J1702-420 A1.47.81.911.0
HESS J1702-420 B6.83.710.5
HESS J1702-420 C8.81.710.5
HESS J1702-420 D0.58.71.410.6
IceCube-170922A0.90.9
Inner Galaxy Survey pos 1-41.03.32.56.8
Inner Galaxy Survey pos 1-50.92.71.95.5
Inner Galaxy Survey pos 1-60.52.32.8
Inner Galaxy Survey pos 1-70.80.90.42.1
Inner Galaxy Survey pos 1-80.51.12.84.4
Inner Galaxy Survey pos 1-91.72.54.2
Inner Galaxy Survey pos 2-50.53.31.95.6
Inner Galaxy Survey pos 2-60.40.4
Inner Galaxy Survey pos 2-74.71.25.9
Inner Galaxy Survey pos 2-81.42.32.15.9
Inner Galaxy Survey pos 3-54.72.87.5
Inner Galaxy Survey pos 3-62.711.60.915.3
Inner Galaxy Survey pos 3-72.31.49.91.915.5
Inner Galaxy Survey pos 3-82.80.58.12.313.7
Kepler SNR2.83.08.214.0
M 870.99.110.0
MAXI J1535-5710.90.9
OJ 2872.32.3
PG 1553+1130.90.30.513.915.6
PKS 0447-4395.31.28.11.315.9
PKS 0754+1001.47.40.50.910.1
PKS 1510-08912.32.160.62.677.5
PKS 2022-0777.87.8
PKS 2155-3043.735.822.81.463.7
PMN 1936-47194.22.26.4
PMN J1548-22510.38.15.50.514.4
PSR B1259-637.00.57.5
Reticulum II1.97.33.312.5
SMC pos 5-113.813.8
SMC pos 5-20.50.5
SMC pos 5-30.50.5
SMC pos 5-40.50.5
SMC pos 5-50.40.4
TXS 0506+0560.410.30.911.6
TXS 0506+0562.82.8
ToO MAXI J1535-5714.72.77.4
ToO Mrk 4211.41.4
ToO PKS 2022-0776.36.3
ToO PKS 2155-3043.90.54.4
ToO170130A1.61.6
ToO1709225.73.69.3
ToO170923A0.20.50.7
ToO171020A3.33.3
ToO20170503A1.61.6
ToOFRB1211020.51.31.8
Tuc III1.39.94.015.2
Tucana II3.26.49.6
V2293 Oph1.10.71.8
Vela Pulsar0.20.2
Vela X PSR2.32.3
Vela pulsar27.62.234.40.965.1
Wd1PosA0.53.64.0
Westerlund 10.50.5
Westerlund 1 Pointing A7.90.526.43.037.8
Westerlund 1 Pointing B13.52.322.55.443.7
Westerlund I0.30.3
custom2.12.1
‘ToO170926’1.91.9
TOTAL195.710.26.9620.2143.34.7981.0
From 59931737dee8b27c5e36655a538c87aee9c2c440 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Thu, 10 Jan 2019 09:34:03 +0100 Subject: [PATCH 2/7] Remove serializer, fixes #887 (#913) --- ctapipe/io/serializer.py | 336 --------------------------- ctapipe/io/tests/test_serializer.py | 155 ------------ docs/development/code-guidelines.rst | 6 - docs/io/index.rst | 10 - examples/obsolete/serialization.py | 59 ----- 5 files changed, 566 deletions(-) delete mode 100644 ctapipe/io/serializer.py delete mode 100644 ctapipe/io/tests/test_serializer.py delete mode 100755 examples/obsolete/serialization.py diff --git a/ctapipe/io/serializer.py b/ctapipe/io/serializer.py deleted file mode 100644 index e16e2415c58..00000000000 --- a/ctapipe/io/serializer.py +++ /dev/null @@ -1,336 +0,0 @@ -""" -Serialize ctapipe containers to file -""" - -from abc import ABC, abstractmethod -from gzip import open as gzip_open -from pickle import dump - -import numpy as np -from astropy import log -from astropy.table import Table, Column -from traitlets import Unicode - -from ctapipe.core import Container - -__all__ = ['Serializer'] - - -class Serializer: - """ - Serializes ctapipe.core.Component, write it to a file thanks - to its Writer object - For some formats (i.e. pickle +gzip), read serialized components from - a file - - Examples - -------- - >>> writer = Serializer(filename='output.pickle', format='pickle', mode='w') - >>> for container in input_containers: - ... writer.add_container(container.r0) - >>> writer.close() - - or using the context manager syntax - >>> with Serializer(filename='output.fits', format='fits', mode='w') as writer: - >>> for container in input_containers: - ... writer.add_container(container.r0) - """ - - def __init__(self, filename, format='fits', mode='x'): - - """ - Parameters - ---------- - filename: str - full path name for i/o file - format: str ('fits', 'img', 'pickle') - mode: str ('write', 'read') - : use this serializer as writer or reader - mode: str - 'w' open for writing, truncating the file first - 'x' open for exclusive creation, failing if the file already exists - 'a' open for writing, appending to the end of the file if it exists - Raises - ------ - NotImplementedError: when format is not implemented - ValueError: when mode is not correct - """ - self.filename = filename - self.format = format - self._stat = None # TODO collect statistics about serialized contents - if mode not in ('x', 'w', 'a'): - raise ValueError('{} is not a valid write mode. Use x, w or a'. - format(mode)) - self._writer = None - - if self.format == 'fits': - self._writer = TableWriter(outfile=filename, mode=mode, - format=format) - elif self.format == 'pickle': - self._writer = GZipPickleWriter(outfile=filename, mode=mode) - elif self.format == 'img': - raise NotImplementedError('img serializer format is' - ' not yet implemented') - else: - raise ValueError('You can serialize only on pickle, fits or img') - - def __enter__(self): - return self - - def __exit__(self, exc_type, exc_value, traceback): - """ - Exit the runtime context related to this object. - The parameters describe the exception that caused the context to be - exited. If the context was exited without an exception, - all three arguments will be None. - If an exception is supplied, and the method wishes to suppress - the exception (i.e., prevent it from being propagated), - it should return a true value. Otherwise, the exception will be - processed normally upon exit from this method. - """ - self.close() - - def add_container(self, container): - """ - Add a container to serializer - """ - self._writer.add_container(container) - - - def close(self): - """ - Write data to disk - """ - self._writer.close() - - -class Writer(ABC): - - def __init__(self, filename): - self.outfile = filename - - @abstractmethod - def add_container(self, container): - pass - - @abstractmethod - def close(self): - pass - - -class GZipPickleWriter(Writer): - """ - Serializes list of ctapipe.core.Components. - Write Component to file - """ - - def __init__(self, outfile, mode='x'): - """ - Parameters - ---------- - outfile: Unicode - full path output file name - mode: str - 'w' open for writing, truncating the file first - 'x' open for exclusive creation, failing if the file already exists - 'a' open for writing, appending to the end of the file if it exists - Raises - ------ - FileNotFoundError: When the file cannot be opened - FileExistsError: when infile exist and mode is x - """ - super().__init__(outfile) - mode += 'b' - try: - self.file_object = gzip_open(outfile, mode) - except FileExistsError: - raise FileExistsError('file exists: {} and mode is {}'. - format(outfile, mode)) - - def close(self): - """ - close opened file - Returns - ------- - """ - self.file_object.close() - - def add_container(self, container): - """ - Add a container to be serialized - - Raises - ------ - TypeError: When container is not type of container - """ - if not isinstance(container, Container): - raise TypeError('Can write only Containers') - dump(container, self.file_object) - - -# FITS Implementation - -not_writeable_fields = ('tel', 'tels_with_data', 'calibration_parameters', - 'pedestal_subtracted_adc', 'integration_window') - - -def is_writeable(key, out_format='fits'): - """ - check if a key is writable - - Parameters - ---------- - key: str - out_format: 'fits' or ´ pickle' - according to out_format a same key can be writable or not - - Returns - ------- - True if key is writable according to the out_format - - Raises - ------ - NameError: When out_format is not know - """ - if out_format is 'fits': - return not (key in not_writeable_fields) - elif out_format is 'pickle': - return True - else: - raise NameError('{} not implemented'.format(out_format)) - - -def writeable_items(container): - """ - # Strip off what we cannot write - Parameters - ---------- - container: ctapipe.core.Container - - Returns - ------- - a dictionary with writable values only - """ - - d = dict(container.items()) - for k in not_writeable_fields: - log.debug("Cannot write column {0}".format(k)) - d.pop(k, None) - return d - - -def to_table(container): - """ - Convert a `ctapipe.core.Container` to an `astropy.Table` with one row - - Parameters - ---------- - container: ctapipe.core.Container - - Returns - ------- - Table: astropy.Table - """ - names = list() - columns = list() - for k, v in writeable_items(container).items(): - - v_arr = np.array(v) - v_arr = v_arr.reshape((1,) + v_arr.shape) - log.debug("Creating column for item '{0}' of shape {1}". - format(k, v_arr.shape)) - names.append(k) - columns.append(Column(v_arr)) - - return Table(data=columns, # dtypes are inferred by columns - names=names, - meta=container.meta) - - -class TableWriter(Writer): - """ - Fits table writer - """ - - def __init__(self, outfile, format='fits', mode='w'): - """ - Parameters - ---------- - outfile: str - output file name - format: str - 'fits' or 'img' - mode: str - 'w' open for writing, truncating the file first - 'x' open for exclusive creation, failing if the file already exists - Raises - ------ - NotImplementedError: when mode is correct but not yet implemented - ValueError: when mode is not correct - """ - super().__init__(outfile) - self.table = Table() - self._created_table = False - self.format = format - self.outfile = outfile - if mode == 'w': - self.overwrite = True - elif mode == 'x': - self.overwrite = False - elif mode == 'a': - raise NotImplementedError('a is a valid write mode,' - ' but not yet implemented') - else: - raise ValueError('{} is not a valid write mode. Use x, w or a'. - format(mode)) - - def _setup_table(self, container): - """ - Create Fits table and HDU - - Parameters - ---------- - container: ctapipe.core.Container - """ - # Create Table from Container - self.table = to_table(container) - - # Write HDU name - if self.format == "fits": - self.table.meta["EXTNAME"] = type(container).__name__ - self._created_table = True - - def add_container(self, container): - """ - Add a container as a table row - Parameters - ---------- - container: ctapipe.core.Container - - Raises - ------ - TypeError: When add another type than Container - """ - if not isinstance(container, Container): - raise TypeError("Can write only Containers") - - if not self._created_table: - self._setup_table(container) - else: - self.table.add_row(writeable_items(container)) - - def close(self, **kwargs): - """ - Write Fits table to file - Parameters - ---------- - kwargs to be passed to `astropy.Table.write method` - - Returns - ------- - Fits Table - """ - # Write table using astropy.table write method - self.table.write(output=self.outfile, format=self.format, - overwrite=self.overwrite, **kwargs) - return self.table diff --git a/ctapipe/io/tests/test_serializer.py b/ctapipe/io/tests/test_serializer.py deleted file mode 100644 index 18d35f14c24..00000000000 --- a/ctapipe/io/tests/test_serializer.py +++ /dev/null @@ -1,155 +0,0 @@ -from copy import deepcopy -from os import remove - -import pytest -from astropy.io import fits - -from ctapipe.io import event_source -from ctapipe.io.serializer import Serializer -from ctapipe.io.sources import PickleSource -from ctapipe.utils import get_dataset_path - - -def compare(read_container, source_container): - # test if 4th adc value of telescope 17 HI_GAIN are equals - return (read_container.r0.tel[17].waveform[0][2][4] == - source_container.r0.tel[17].waveform[0][2][4]) - - -def generate_input_containers(): - # Get event from hessio file, append them into input_containers - input_filename = get_dataset_path("gamma_test.simtel.gz") - with event_source(input_filename, max_events=3) as source: - input_containers = [deepcopy(event) for event in source] - return input_containers - - -# Setup -input_containers = generate_input_containers() - - -@pytest.fixture(scope='session') -def binary_filename(tmpdir_factory): - return str(tmpdir_factory.mktemp('data') - .join('pickle_data.pickle.gz')) - - -@pytest.fixture(scope='session') -def fits_file_name(tmpdir_factory): - return str(tmpdir_factory.mktemp('data').join('output.fits')) - - -def test_pickle_serializer(binary_filename): - serial = Serializer(filename=binary_filename, format='pickle', mode='w') - # append all input file events in input_containers list and pickle serializer - for event in input_containers: - serial.add_container(event) - serial.close() - - # read Containers from pickle serializer - reader = PickleSource(filename=binary_filename) - # file read_containers from serializer generator - read_containers = [] - for container in reader: - read_containers.append(container) - # test if number of read Container correspond to input - assert len(read_containers) is len(input_containers) - # test if 4th adc value of telescope 17 HI_GAIN are equals - assert compare(input_containers[2], read_containers[2]) - reader.close() - remove(binary_filename) - - -# Test pickle reader/writer with statement -def test_pickle_with_statement(binary_filename): - with Serializer(filename=binary_filename, format='pickle', mode='w') as \ - containers_writer: - for container in input_containers: - containers_writer.add_container(container) - containers_writer.close() - - read_containers = [] - with PickleSource(filename=binary_filename) as reader: - for container in reader: - read_containers.append(container) - # test if number of read Container correspond to input - assert len(read_containers) is len(input_containers) - # test if 4th adc value of telescope 17 HI_GAIN are equals - assert compare(input_containers[2], read_containers[2]) - remove(binary_filename) - - -# Test pickle reader iterator -def test_pickle_iterator(binary_filename): - serial = Serializer(filename=binary_filename, format='pickle', - mode='w') - # append all events in input_containers list and pickle serializer - for event in input_containers: - serial.add_container(event) - serial.close() - - read_containers = [] - reader = PickleSource(filename=binary_filename) - for container in reader: - read_containers.append(container) - # test if number of read Container correspond to input - assert len(read_containers) is len(input_containers) - # test if 4th adc value of telescope 17 HI_GAIN are equals - assert compare(input_containers[2], read_containers[2]) - reader.close() - remove(binary_filename) - - - - - -def test_fits_dl0(fits_file_name): - serial = Serializer(filename=fits_file_name, format='fits', mode='w') - for container in input_containers: - serial.add_container(container.dl0) - serial.close() - hdu = fits.open(fits_file_name)[1] - assert hdu.data["event_id"][0] == 408 - assert hdu.data["event_id"][1] == 409 - assert hdu.data["event_id"][2] == 803 - assert hdu.data["obs_id"][2] == 31964 - remove(fits_file_name) - - -def test_exclusive_mode(fits_file_name): - serial = Serializer(filename=fits_file_name, format='fits', mode='w') - for container in input_containers: - serial.add_container(container.dl0) - serial.close() - # Try to write to fits_file_name in exclusive mode - with pytest.raises(OSError): - serial = Serializer(filename=fits_file_name, format='fits', mode='x') - serial.add_container(input_containers[2].dl0) - serial.close() - remove(fits_file_name) - -""" -def test_fits_dl1(): - input_test_file = get_datasets_path('example_container.pickle.gz') - with gzip_open(input_test_file, 'rb') as f: - data = load(f) - t38 = data[0].dl1.tel[38] - serial = Serializer('output.fits', 'fits', overwrite=True) - serial.add_container(t38) - serial.write() - # t11_1 = data[1].dl1.tel[11] - # S_cal.write(t11_1) # This will not work because shape of data is different from tel to tel. -""" - - -def test_fits_context_manager(fits_file_name): - with Serializer(filename=fits_file_name, format='fits', mode='w') as writer: - for container in input_containers: - writer.add_container(container.dl0) - - hdulist = fits.open(fits_file_name) - assert hdulist[1].data["event_id"][0] == 408 - remove(fits_file_name) - - -# TODO test FITSSource class diff --git a/docs/development/code-guidelines.rst b/docs/development/code-guidelines.rst index ca08c29c5a6..915313893c8 100644 --- a/docs/development/code-guidelines.rst +++ b/docs/development/code-guidelines.rst @@ -306,7 +306,6 @@ help when writing algorithms: source = calibrated_event_source(filename) ImageMangler mangler(geom.pix_x, geom.pix_y, "transformtable.fits") - Serializer serializer = ... # simple loop over events, calling each algorithm and directly #passing data @@ -317,11 +316,6 @@ help when writing algorithms: mangled_image = mangler.mangle(image) image_parameters = parameterize_image(mangled_image) - # here you may here pack your output values into a Container if - # they are not already in one. We assume here that mangled_image - # and image_parameters are already Container subclasses - - serializer.write([mangled_image, image_parameters]) * When your algorithm test code (as above) works well and you are happy with the results, you can do two things: diff --git a/docs/io/index.rst b/docs/io/index.rst index a7e0e07f7d7..3b6d79391a6 100644 --- a/docs/io/index.rst +++ b/docs/io/index.rst @@ -160,10 +160,6 @@ data: Serialization of Containers: ============================ -The `serializer` module provide support for storing -`ctapipe.io.Container` classes in output files (for example FITS -tables or pickle files) - The `ctapipe.io.TableWriter` and `ctapipe.io.TableReader` base classes provide an interface to implement subclasses that write/read Containers to/from table-like data files. Currently the only implementation is for writing @@ -188,9 +184,3 @@ Reference/API .. automodapi:: ctapipe.io.containers :no-inheritance-diagram: - ------------------------------- - -.. automodapi:: ctapipe.io.serializer - :no-inheritance-diagram: - diff --git a/examples/obsolete/serialization.py b/examples/obsolete/serialization.py deleted file mode 100755 index dc2d1d1caee..00000000000 --- a/examples/obsolete/serialization.py +++ /dev/null @@ -1,59 +0,0 @@ -#!/usr/bin/env python3 - -from ctapipe.io.serializer import Serializer -from astropy import log -import pickle -import gzip - -log.setLevel('DEBUG') - - -def write_dl0_example(filename, data): - S = Serializer(filename, mode='w') - - # Create table - for container in data: - S.write(container.dl0) - - # Print table - print(S._writer.table) - - # Save table to disk - S.save() - return S - - -def write_dl1_tel_example(filename, data): - - t38 = data[0].dl1.tel[38] - - S_cal = Serializer(filename, mode='w') - S_cal.write(t38) - - print(S_cal._writer.table) - - # t11_1 = data[1].dl1.tel[11] - # S_cal.write(t11_1) - # This will not work because shape of data is different from tel to tel. - - S_cal.save() - return S_cal - - -def context_manager_example(filename, data): - with Serializer(filename, mode='w') as writer: - for container in data: - print(container.dl0) - writer.write(container.dl0) - print(writer._writer.table) - return 0 - - -if __name__ == "__main__": - - with gzip.open('example_container.pickle.gz', 'rb') as f: - data = pickle.load(f) - - S = write_dl0_example("output.fits", data) - S_cal = write_dl1_tel_example("cal_output.fits", data) - S_context = context_manager_example("output_context.fits", data) From 097f0e04787d8c2e133075b74995f5e4f1ab8709 Mon Sep 17 00:00:00 2001 From: Jason Watson Date: Fri, 11 Jan 2019 17:47:49 +0100 Subject: [PATCH 3/7] Charge Resolution Update (#827) * Updated charge resolution handling - Renamed ctapipe/analysis/camera/chargeresolution.py to ctapipe/analysis/camera/charge_resolution.py - Updated ChargeResolutionCalculator to be much simpler. Now creates a list of pandas dataframes, which are combined once the list is above 1GB in memory. Also now allows non-integer true charge values (useful for lab measurements). Charge resolution is then stored in a HDF5 file via a pandas HDFStore, a much simpler interface than the pre-existing storage format for charge resolutions. - Created ctapipe/plotting/charge_resolution.py to handle the plotting of the charge resolution files. - Updated tools to use the new classes * Updated tests for charge resolution classes * Removed left-over comments * Deleted plot_charge_resolution_variation_hist * Removed ctapipe-chargeres-hist from setup.py * Updated description for plot_charge_resolution * Lack of punctuation causing test to fail... * Close figure after saving to fix test failing * Corrected the use of tmpdir Convert LocalPath to str * Further tmpdir corrections to str * Further tmpdir corrections to str * Corrected variable name * Address comments * Add provenance --- ctapipe/analysis/camera/__init__.py | 2 +- ctapipe/analysis/camera/charge_resolution.py | 135 +++++++ ctapipe/analysis/camera/chargeresolution.py | 272 ------------- .../camera/tests/test_charge_resolution.py | 100 +++++ .../camera/tests/test_chargeresolution.py | 67 ---- ctapipe/plotting/charge_resolution.py | 367 ++++++++++++++++++ .../plotting/tests/test_charge_resolution.py | 124 ++++++ ctapipe/tools/extract_charge_resolution.py | 127 +++--- ctapipe/tools/plot_charge_resolution.py | 195 ++-------- .../plot_charge_resolution_variation_hist.py | 146 ------- ctapipe/tools/tests/test_tools.py | 30 ++ setup.py | 2 - 12 files changed, 844 insertions(+), 723 deletions(-) create mode 100644 ctapipe/analysis/camera/charge_resolution.py create mode 100644 ctapipe/analysis/camera/tests/test_charge_resolution.py delete mode 100644 ctapipe/analysis/camera/tests/test_chargeresolution.py create mode 100644 ctapipe/plotting/charge_resolution.py create mode 100644 ctapipe/plotting/tests/test_charge_resolution.py delete mode 100644 ctapipe/tools/plot_charge_resolution_variation_hist.py diff --git a/ctapipe/analysis/camera/__init__.py b/ctapipe/analysis/camera/__init__.py index fe880556385..7a8fcd15b52 100644 --- a/ctapipe/analysis/camera/__init__.py +++ b/ctapipe/analysis/camera/__init__.py @@ -2,4 +2,4 @@ Analysis modules specifically related to camera properties. """ -from .chargeresolution import * +from .charge_resolution import * diff --git a/ctapipe/analysis/camera/charge_resolution.py b/ctapipe/analysis/camera/charge_resolution.py new file mode 100644 index 00000000000..832c161f641 --- /dev/null +++ b/ctapipe/analysis/camera/charge_resolution.py @@ -0,0 +1,135 @@ +import numpy as np +import pandas as pd + +__all__ = ['ChargeResolutionCalculator'] + + +class ChargeResolutionCalculator: + def __init__(self, mc_true=True): + """ + Calculates the charge resolution with an efficient, low-memory, + interative approach, allowing the contribution of data/events + without reading the entire dataset into memory. + + Utilises Pandas DataFrames, and makes no assumptions on the order of + the data, and does not require the true charge to be integer (as may + be the case for lab measurements where an average illumination + is used). + + A list is filled with a dataframe for each contribution, and only + amalgamated into a single dataframe (reducing memory) once the memory + of the list becomes large (or at the end of the filling), + reducing the time required to produce the output. + + Parameters + ---------- + mc_true : bool + Indicate if the "true charge" values are from the sim_telarray + files, and therefore without poisson error. The poisson error will + therefore be included in the charge resolution calculation. + + Attributes + ---------- + self._mc_true : bool + self._df_list : list + self._df : pd.DataFrame + self._n_bytes : int + Monitors the number of bytes being held in memory + """ + self._mc_true = mc_true + self._df_list = [] + self._df = pd.DataFrame() + self._n_bytes = 0 + self._max_bytes = 1E9 + + @staticmethod + def rmse_abs(sum_, n): + return np.sqrt(sum_ / n) + + @staticmethod + def rmse(true, sum_, n): + return ChargeResolutionCalculator.rmse_abs(sum_, n) / np.abs(true) + + @staticmethod + def charge_res_abs(true, sum_, n): + return np.sqrt((sum_ / n) + true) + + @staticmethod + def charge_res(true, sum_, n): + return (ChargeResolutionCalculator.charge_res_abs(true, sum_, n) + / np.abs(true)) + + def add(self, pixel, true, measured): + """ + Contribute additional values to the Charge Resolution + + Parameters + ---------- + pixel : ndarray + 1D array containing the pixel for each entry + true : ndarray + 1D array containing the true charge for each entry + measured : ndarray + 1D array containing the measured charge for each entry + """ + diff2 = np.power(measured - true, 2) + df = pd.DataFrame(dict( + pixel=pixel, + true=true, + sum=diff2, + n=np.uint32(1) + )) + self._df_list.append(df) + self._n_bytes += df.memory_usage(index=True, deep=True).sum() + if self._n_bytes > self._max_bytes: + self._amalgamate() + + def _amalgamate(self): + """ + Concatenate the dataframes inside the list, and sum together + values per pixel and true charge in order to reduce memory use. + """ + self._df = pd.concat([self._df, *self._df_list], ignore_index=True) + self._df = self._df.groupby(['pixel', 'true']).sum().reset_index() + self._n_bytes = 0 + self._df_list = [] + + def finish(self): + """ + Perform the final amalgamation, and calculate the charge resolution + from the resulting sums + + Returns + ------- + df_p : pd.DataFrame + Dataframe containing the charge resolution per pixel + df_c : pd.DataFrame + Dataframe containing the charge resolution for the entire camera + """ + self._amalgamate() + + self._df = self._df.loc[self._df['true'] != 0] + + df_p = self._df.copy() + true = df_p['true'].values + sum_ = df_p['sum'].values + n = df_p['n'].values + if self._mc_true: + df_p['charge_resolution'] = self.charge_res(true, sum_, n) + df_p['charge_resolution_abs'] = self.charge_res_abs(true, sum_, n) + else: + df_p['charge_resolution'] = self.rmse(true, sum_, n) + df_p['charge_resolution_abs'] = self.rmse_abs(sum_, n) + df_c = self._df.copy().groupby('true').sum().reset_index() + df_c = df_c.drop(columns='pixel') + true = df_c['true'].values + sum_ = df_c['sum'].values + n = df_c['n'].values + if self._mc_true: + df_c['charge_resolution'] = self.charge_res(true, sum_, n) + df_c['charge_resolution_abs'] = self.charge_res_abs(true, sum_, n) + else: + df_c['charge_resolution'] = self.rmse(true, sum_, n) + df_c['charge_resolution_abs'] = self.rmse_abs(sum_, n) + + return df_p, df_c diff --git a/ctapipe/analysis/camera/chargeresolution.py b/ctapipe/analysis/camera/chargeresolution.py index a901154d49f..e69de29bb2d 100644 --- a/ctapipe/analysis/camera/chargeresolution.py +++ b/ctapipe/analysis/camera/chargeresolution.py @@ -1,272 +0,0 @@ -from math import log10, sqrt -from os import makedirs -from os.path import dirname, exists - -import numpy as np -from scipy.stats import binned_statistic as bs -from tables import open_file - -from ctapipe.core import Component -from ctapipe.core.traits import Int, Bool - -__all__ = ['ChargeResolutionCalculator'] - - -class ChargeResolutionCalculator(Component): - """ - Class to handle the calculation of Charge Resolution. - - Attributes - ---------- - max_pe : int - Maximum pe to calculate the charge resolution up to. - sum_dict : dict - Dictionary to store the running sum for each true charge. - n_dict : dict - Dictionary to store the running number for each true charge. - variation_hist_nbins : float - Number of bins for the variation histogram. - variation_hist_range : list - X and Y range for the variation histogram. - variation_hist : `np.histogram2d` - variation_xedges : ndarray - Edges of the X bins for the variation histogram. - variation_yedges : ndarray - Edges of the Y bins for the variation histogram. - """ - max_pe = Int(2000, help='Maximum pe to calculate the charge resolution ' - 'up to').tag(config=True) - binning = Int(60, allow_none=True, - help='Number of bins for the Charge Resolution. If None, ' - 'no binning is performed.').tag(config=True) - log_bins = Bool(True, help='Bin the x axis linearly instead of ' - 'logarithmic.').tag(config=True) - - def __init__(self, config=None, tool=None, **kwargs): - """ - Calculator of charge resolution. - - Parameters - ---------- - config : traitlets.loader.Config - Configuration specified by config file or cmdline arguments. - Used to set traitlet values. - Set to None if no configuration to pass. - tool : ctapipe.core.Tool - Tool executable that is calling this component. - Passes the correct logger to the component. - Set to None if no Tool to pass. - kwargs - """ - super().__init__(config=config, parent=tool, **kwargs) - - self.sum_array = np.zeros(self.max_pe) - self.n_array = np.zeros(self.max_pe) - self.sum_dict = {} - self.n_dict = {} - - self.variation_hist_nbins = int(log10(self.max_pe) * 50) - self.variation_hist_range = [[log10(1), log10(self.max_pe)], - [log10(1), log10(self.max_pe)]] - h, xedges, yedges = np.histogram2d([np.nan], [np.nan], - bins=self.variation_hist_nbins, - range=self.variation_hist_range) - self.variation_hist = h - self.variation_xedges = xedges - self.variation_yedges = yedges - - self.storage_arrays = ['max_pe', 'sum_array', 'n_array', - 'variation_hist_nbins', 'variation_hist_range', - 'variation_hist', 'variation_xedges', - 'variation_yedges'] - - def add_charges(self, true_charge, measured_charge): - """ - Fill the class parameters with a numpy array of true charge and - measured (calibrated) charge from an event. The two arrays must be the - same size. - - Parameters - ---------- - true_charge : ndarray - Array of true (MC) charge. - Obtained from event.mc.tel[telid].image[channel] - measured_charge : ndarray - Array of measured (dl1 calibrated) charge. - Obtained from event.mc.tel[tel_id].photo_electron_image - """ - above_0 = (measured_charge > 0) & (true_charge > 0) - x = true_charge[above_0] - y = measured_charge[above_0] - h, _, _ = np.histogram2d(np.log10(y), - np.log10(x), - bins=self.variation_hist_nbins, - range=self.variation_hist_range) - self.variation_hist += h - - in_range = (true_charge > 0) & (true_charge <= self.max_pe) - true_q = true_charge[in_range] - measured_q = measured_charge[in_range] - np.add.at(self.sum_array, true_q - 1, np.power(measured_q - true_q, 2)) - np.add.at(self.n_array, true_q - 1, 1) - - def get_charge_resolution(self): - """ - Calculate and obtain the charge resolution graph arrays. - - Returns - ------- - true_charge : ndarray - The X axis true charges. - chargeres : ndarray - The Y axis charge resolution values. - chargeres_error : ndarray - The error on the charge resolution. - scaled_chargeres : ndarray - The Y axis charge resolution divided by the Goal. - scaled_chargeres_error : ndarray - The error on the charge resolution divided by the Goal. - """ - self.log.debug('[chargeres] Calculating charge resolution') - - n_1 = self.n_array > 0 - n = self.n_array[n_1] - true = (np.arange(self.max_pe) + 1)[n_1] - sum_ = self.sum_array[n_1] - - res = np.sqrt((sum_ / n) + true) / true - res_error = res * (1 / np.sqrt(2 * n)) - - scale = self.goal(true) - scaled_res = res / scale - scaled_res_error = res_error / scale - - if self.binning is not None: - x = true - if self.log_bins: - x = np.log10(true) - - def binning(array): - return bs(x, array, 'mean', bins=self.binning) - - def sum_errors(array): - return np.sqrt(np.sum(np.power(array, 2))) / array.size - - def bin_errors(array): - return bs(x, array, sum_errors, bins=self.binning) - - true, _, _ = binning(true) - res, _, _ = binning(res) - res_error, _, _ = bin_errors(res_error) - scaled_res, _, _ = binning(scaled_res) - scaled_res_error, _, _ = bin_errors(scaled_res_error) - - return true, res, res_error, scaled_res, scaled_res_error - - @staticmethod - def limit_curves(npe, n_nsb, n_add, enf, sigma2): - """ - Equation for calculating the Goal and Requirement curves, as defined - in SCI-MC/121113. - https://portal.cta-observatory.org/recordscentre/Records/SCI/ - SCI-MC/measurment_errors_system_performance_1YQCBC.pdf - - Parameters - ---------- - npe : ndarray - Number of photoeletrons (variable). - n_nsb : float - Number of NSB photons. - n_add : float - Number of photoelectrons from additional noise sources. - enf : float - Excess noise factor. - sigma2 : float - Percentage ofmultiplicative errors. - """ - return (np.sqrt((n_nsb + n_add) + np.power(enf, 2) * npe + - np.power(sigma2 * npe, 2)) / npe).astype(float) - - @staticmethod - def requirement(npe): - """ - CTA requirement curve. - - Parameters - ---------- - npe : ndarray - Number of photoeletrons (variable). - """ - n_nsb = sqrt(4.0 + 3.0) - n_add = 0 - enf = 1.2 - sigma2 = 0.1 - defined_npe = 1000 - - # If npe is not an array, temporarily convert it to one - npe = np.array([npe]) - - lc = ChargeResolutionCalculator.limit_curves - requirement = lc(npe, n_nsb, n_add, enf, sigma2) - requirement[npe > defined_npe] = np.nan - - return requirement[0] - - @staticmethod - def goal(npe): - """ - CTA goal curve. - - Parameters - ---------- - npe : ndarray - Number of photoeletrons (variable). - """ - n_nsb = 2 - n_add = 0 - enf = 1.1152 - sigma2 = 0.05 - defined_npe = 2000 - - # If npe is not an array, temporarily convert it to one - npe = np.array([npe]) - - lc = ChargeResolutionCalculator.limit_curves - goal = lc(npe, n_nsb, n_add, enf, sigma2) - goal[npe > defined_npe] = np.nan - - return goal[0] - - @staticmethod - def poisson(npe): - """ - Poisson limit curve. - - Parameters - ---------- - npe : ndarray - Number of photoeletrons (variable). - """ - # If npe is not an array, temporarily convert it to one - npe = np.array([npe]) - poisson = np.sqrt(npe) / npe - - return poisson[0] - - def save(self, path): - output_dir = dirname(path) - if not exists(output_dir): - self.log.info("[output] Creating directory: {}".format(output_dir)) - makedirs(output_dir) - self.log.info("Saving Charge Resolution file: {}".format(path)) - - with open_file(path, mode="w", title="ChargeResolutionFile") as f: - group = f.create_group("/", 'ChargeResolution', '') - for arr in self.storage_arrays: - f.create_array(group, arr, getattr(self, arr), arr) - - def load(self, path): - self.log.info("Loading Charge Resolution file: {}".format(path)) - with open_file(path, mode="r") as f: - for arr in self.storage_arrays: - setattr(self, arr, f.get_node("/ChargeResolution", arr).read()) diff --git a/ctapipe/analysis/camera/tests/test_charge_resolution.py b/ctapipe/analysis/camera/tests/test_charge_resolution.py new file mode 100644 index 00000000000..decc5bc8c56 --- /dev/null +++ b/ctapipe/analysis/camera/tests/test_charge_resolution.py @@ -0,0 +1,100 @@ +from ctapipe.analysis.camera.charge_resolution import \ + ChargeResolutionCalculator +import numpy as np +from numpy.testing import assert_almost_equal + + +def test_add(): + chargeres = ChargeResolutionCalculator() + true_charge = np.arange(100) + measured_charge = np.arange(100) + chargeres.add(0, true_charge, measured_charge) + assert len(chargeres._df_list) == 1 + assert chargeres._df_list[0].index.size == 100 + + +def test_amalgamate(): + chargeres = ChargeResolutionCalculator() + true_charge = np.arange(100) + measured_charge = np.arange(100) + chargeres.add(0, true_charge, measured_charge) + chargeres.add(0, true_charge, measured_charge) + chargeres.add(1, true_charge, measured_charge) + assert len(chargeres._df_list) == 3 + chargeres._amalgamate() + assert len(chargeres._df_list) == 0 + assert chargeres._df.index.size == 200 + + +def test_memory_limit(): + chargeres = ChargeResolutionCalculator() + chargeres._max_bytes = 1 + true_charge = np.arange(100) + measured_charge = np.arange(100) + chargeres.add(0, true_charge, measured_charge) + assert len(chargeres._df_list) == 0 + chargeres.add(0, true_charge, measured_charge) + assert len(chargeres._df_list) == 0 + + +def test_finish(): + chargeres = ChargeResolutionCalculator() + true_charge = np.arange(100) + measured_charge = np.arange(100) + chargeres.add(0, true_charge, measured_charge) + chargeres.add(0, true_charge, measured_charge) + df_p, df_c = chargeres.finish() + assert np.array_equal(df_p['charge_resolution'].values, + df_c['charge_resolution'].values) + chargeres.add(1, true_charge, measured_charge) + df_p, df_c = chargeres.finish() + assert not np.array_equal(df_p['charge_resolution'].values, + df_c['charge_resolution'].values) + + +def test_calculation(): + chargeres = ChargeResolutionCalculator() + measured = np.array([3.5, 2.7]) + true = 3 + n = measured.size + + sum_ = np.sum(np.power(measured - true, 2)) + assert_almost_equal(sum_, 0.34, 3) + assert_almost_equal(chargeres.rmse_abs(sum_, n), 0.412, 3) + assert_almost_equal(chargeres.rmse(true, sum_, n), 0.137, 3) + assert_almost_equal(chargeres.charge_res_abs(true, sum_, n), 1.780, 3) + assert_almost_equal(chargeres.charge_res(true, sum_, n), 0.593, 3) + + assert chargeres.rmse_abs(sum_, n) == chargeres.rmse(true, sum_, n) * true + assert (chargeres.charge_res_abs(true, sum_, n) == + chargeres.charge_res(true, sum_, n) * true) + + +def test_result(): + chargeres = ChargeResolutionCalculator(mc_true=False) + measured = np.array([3.5, 2.7]) + true = 3 + n = measured.size + sum_ = np.sum(np.power(measured - true, 2)) + + chargeres.add(0, true, measured) + df_p, df_c = chargeres.finish() + assert (df_p['charge_resolution'].values[0] == + chargeres.rmse(true, sum_, n)) + assert (df_p['charge_resolution_abs'].values[0] == + chargeres.rmse_abs(sum_, n)) + + +def test_result_mc_true(): + chargeres = ChargeResolutionCalculator() + measured = np.array([3.5, 2.7]) + true = 3 + n = measured.size + sum_ = np.sum(np.power(measured - true, 2)) + + chargeres.add(0, true, measured) + df_p, df_c = chargeres.finish() + assert (df_p['charge_resolution'].values[0] == + chargeres.charge_res(true, sum_, n)) + assert (df_p['charge_resolution_abs'].values[0] == + chargeres.charge_res_abs(true, sum_, n)) diff --git a/ctapipe/analysis/camera/tests/test_chargeresolution.py b/ctapipe/analysis/camera/tests/test_chargeresolution.py deleted file mode 100644 index 8a733edcd03..00000000000 --- a/ctapipe/analysis/camera/tests/test_chargeresolution.py +++ /dev/null @@ -1,67 +0,0 @@ -from ctapipe.analysis.camera.chargeresolution import ChargeResolutionCalculator -import numpy as np - - -def test_add_charges(): - chargeres = ChargeResolutionCalculator() - true_charge = np.arange(100) - measured_charge = np.arange(100) - chargeres.add_charges(true_charge, measured_charge) - variation_hist = chargeres.variation_hist - - assert(variation_hist[0][0] == 1) - - -# def test_add_source(): - # Cannot currently test - no test file has true charge - - -def test_get_charge_resolution(): - chargeres = ChargeResolutionCalculator( binning=None) - true_charge = np.arange(100) - measured_charge = np.arange(100) - chargeres.add_charges(true_charge, measured_charge) - true_charge, chargeres, chargeres_error, \ - scaled_chargeres, scaled_chargeres_error = \ - chargeres.get_charge_resolution() - - assert(true_charge[0] == 1) - assert(chargeres[0] == 1) - assert(round(chargeres_error[0], 7) == 0.7071068) - assert(round(scaled_chargeres[0], 7) == 0.5550272) - assert(round(scaled_chargeres_error[0], 7) == 0.3924635) - - -def test_get_binned_charge_resolution(): - chargeres = ChargeResolutionCalculator() - true_charge = np.arange(100) - measured_charge = np.arange(100) - chargeres.add_charges(true_charge, measured_charge) - true_charge, chargeres, chargeres_error, \ - scaled_chargeres, scaled_chargeres_error = \ - chargeres.get_charge_resolution() - - assert(true_charge[0] == 1) - assert(chargeres[0] == 1) - assert(round(chargeres_error[0], 7) == 0.7071068) - assert(round(scaled_chargeres[0], 7) == 0.5550272) - assert(round(scaled_chargeres_error[0], 7) == 0.3924635) - - -def test_limit_curves(): - value = ChargeResolutionCalculator.limit_curves(1, 2, 3, 4, 5) - assert(round(value, 7) == 6.78233) - - -def test_requirement(): - value = ChargeResolutionCalculator.requirement(1) - value2 = ChargeResolutionCalculator.requirement(np.arange(1, 10))[-1] - assert(round(value, 7) == 2.0237963) - assert(round(value2, 7) == 0.4501817) - - -def test_goal(): - value = ChargeResolutionCalculator.goal(1) - value2 = ChargeResolutionCalculator.goal(np.arange(1, 10))[-1] - assert(round(value, 7) == 1.8017134) - assert(round(value2, 7) == 0.4066657) diff --git a/ctapipe/plotting/charge_resolution.py b/ctapipe/plotting/charge_resolution.py new file mode 100644 index 00000000000..1adc4ab4359 --- /dev/null +++ b/ctapipe/plotting/charge_resolution.py @@ -0,0 +1,367 @@ +import numpy as np +import pandas as pd +from matplotlib.ticker import FuncFormatter +from matplotlib import pyplot as plt +import os +from ctapipe.core import Component, Provenance +from traitlets import Unicode, Int + + +def sum_errors(array): + """ + Simple sum of squares to combine errors + + Parameters + ---------- + array : ndarray + + Returns + ------- + float + """ + return np.sqrt(np.sum(np.power(array, 2)) / array.size) + + +def bin_dataframe(df, n_bins): + """ + Assign a "bin" column to the dataframe to indicate which bin the true + charges belong to. + + Bins are assigned in log space. + + Parameters + ---------- + df : pd.DataFrame + n_bins : int + Number of bins to allow in range + + Returns + ------- + pd.DataFrame + """ + true = df['true'].values + min_ = true.min() + max_ = true.max() + bins = np.geomspace(min_, max_, n_bins) + bins = np.append(bins, 10**(np.log10(bins[-1]) + + np.diff(np.log10(bins))[0])) + df['bin'] = np.digitize(true, bins, right=True) - 1 + + return df + + +class ChargeResolutionPlotter(Component): + + output_path = Unicode( + '', help='Output path to save the plot.' + ).tag(config=True) + + n_bins = Int( + 40, + help='Number of bins for collecting true charges and combining ' + 'their resolution' + ).tag(config=True) + + def __init__(self, config=None, tool=None, **kwargs): + """ + Plots the charge resolution HDF5 file obtained from + `ctapipe.analysis.camera.charge_resolution`. + + Also contains the equation from which the charge resolution + requirement is defined, which can be plotted alongside the charge + resolution curve. + + `ctapipe.tools.plot_charge_resolution` demonstrated the use of this + component. + + Parameters + ---------- + config : traitlets.loader.Config + Configuration specified by config file or cmdline arguments. + Used to set traitlet values. + Set to None if no configuration to pass. + tool : ctapipe.core.Tool + Tool executable that is calling this component. + Passes the correct logger to the component. + Set to None if no Tool to pass. + kwargs + """ + super().__init__(config=config, parent=tool, **kwargs) + self._df_pixel = None + self._df_camera = None + + self.fig = plt.figure(figsize=(12, 7.42)) + self.ax = self.fig.add_subplot(111) + + self.ax.set_xlabel("True Charge (p.e.)") + self.ax.set_ylabel( + r"Fractional Charge Resolution $\frac{{\sigma_Q}}{{Q}}$" + ) + + if not self.output_path: + raise ValueError("Output path must be specified") + + def _set_file(self, path): + """ + Reads the charge resolution DataFrames from the file and assigns it to + this class ready for plotting. + + Parameters + ---------- + path : str + Path to the charge resolution HDF5 file + """ + with pd.HDFStore(path, 'r') as store: + self._df_pixel = store['charge_resolution_pixel'] + self._df_camera = store['charge_resolution_camera'] + + def _plot(self, x, y, **kwargs): + """ + Plot the given points onto the figure + + Parameters + ---------- + x : ndarray + y : ndarray + xerr : ndarray + yerr : ndarray + label : str + """ + defaults = dict( + mew=1, capsize=1, elinewidth=0.5, markersize=2, + linewidth=0.5, fmt='.' + ) + kwargs = {**defaults, **kwargs} + (_, caps, _) = self.ax.errorbar(x, y, **kwargs) + for cap in caps: + cap.set_markeredgewidth(0.5) + + def plot_average(self, path, label='', **kwargs): + """ + Plot the average and standard deviation of the charge resolution + across the pixels of the camera. + + Parameters + ---------- + path : str + Path to the charge resolution HDF5 file + label : str + Label for the figure's legend + kwargs + """ + self._set_file(path) + df_binned = bin_dataframe(self._df_pixel, self.n_bins) + agg = {'charge_resolution': ['mean', 'std'], 'true': 'mean'} + df_agg = df_binned.groupby(['bin']).agg(agg) + x = df_agg['true']['mean'].values + y = df_agg['charge_resolution']['mean'].values + yerr = df_agg['charge_resolution']['std'].values + self._plot(x, y, yerr=yerr, label=label, **kwargs) + + def plot_pixel(self, path, pixel, label='', **kwargs): + """ + Plot a single pixel's charge resolution. + + The yerr represents the amount of entries. + + Parameters + ---------- + path : str + Path to the charge resolution HDF5 file + pixel : int + Pixel index to plot + label : str + Label for the figure's legend + kwargs + """ + self._set_file(path) + df_p = self._df_pixel.loc[self._df_pixel['pixel'] == pixel] + df_binned = bin_dataframe(df_p, self.n_bins) + agg = {'charge_resolution': 'mean', 'true': 'mean', 'n': 'sum'} + df_agg = df_binned.groupby(['bin']).agg(agg) + x = df_agg['true'].values + y = df_agg['charge_resolution'].values + yerr = 1 / np.sqrt(df_agg['n'].values) + self._plot(x, y, yerr=yerr, label=label, **kwargs) + + def plot_camera(self, path, label='', **kwargs): + """ + Plot the charge resolution for the entire camera. + + The yerr represents the amount of entries. + + Parameters + ---------- + path : str + Path to the charge resolution HDF5 file + label : str + Label for the figure's legend + kwargs + """ + self._set_file(path) + df_binned = bin_dataframe(self._df_camera, self.n_bins) + agg = {'charge_resolution': 'mean', 'true': 'mean', 'n': 'sum'} + df_agg = df_binned.groupby(['bin']).agg(agg) + x = df_agg['true'].values + y = df_agg['charge_resolution'].values + yerr = 1 / np.sqrt(df_agg['n'].values) + self._plot(x, y, yerr=yerr, label=label, **kwargs) + + def _finish(self): + """ + Perform the finishing touches to the figure before saving. + """ + self.ax.set_xscale('log') + self.ax.get_xaxis().set_major_formatter( + FuncFormatter(lambda x, _: '{:g}'.format(x))) + self.ax.set_yscale('log') + self.ax.get_yaxis().set_major_formatter( + FuncFormatter(lambda y, _: '{:g}'.format(y))) + + def save(self): + """ + Save the figure to the path defined by the `output_path` trait + """ + self._finish() + + output_dir = os.path.dirname(self.output_path) + if not os.path.exists(output_dir): + print("Creating directory: {}".format(output_dir)) + os.makedirs(output_dir) + + self.fig.savefig(self.output_path, bbox_inches='tight') + print("Figure saved to: {}".format(self.output_path)) + Provenance().add_output_file(self.output_path) + + plt.close(self.fig) + + @staticmethod + def limit_curves(q, nsb, t_w, n_e, sigma_g, enf): + """ + Equation for calculating the Goal and Requirement curves, as defined + in SCI-MC/121113. + https://portal.cta-observatory.org/recordscentre/Records/SCI/ + SCI-MC/measurment_errors_system_performance_1YQCBC.pdf + + Parameters + ---------- + q : ndarray + Number of photoeletrons (variable). + nsb : float + Number of NSB photons. + t_w : float + Effective signal readout window size. + n_e : float + Electronic noise + sigma_g : float + Multiplicative errors of the gain. + enf : float + Excess noise factor. + """ + sigma_0 = np.sqrt(nsb * t_w + n_e**2) + sigma_enf = 1 + enf + sigma_q = np.sqrt(sigma_0**2 + sigma_enf**2 * q + sigma_g**2 * q**2) + return sigma_q / q + + @staticmethod + def requirement(q): + """ + CTA requirement curve. + + Parameters + ---------- + q : ndarray + Number of photoeletrons + """ + nsb = 0.125 + t_w = 15 + n_e = 0.87 + sigma_g = 0.1 + enf = 0.2 + defined_npe = 1000 + lc = __class__.limit_curves + requirement = lc(q, nsb, t_w, n_e, sigma_g, enf) + requirement[q > defined_npe] = np.nan + + return requirement + + @staticmethod + def poisson(q): + """ + Poisson limit curve. + + Parameters + ---------- + q : ndarray + Number of photoeletrons + """ + poisson = np.sqrt(q) / q + return poisson + + def plot_requirement(self, q): + """ + Plot the CTA requirement curve onto the figure. + + Parameters + ---------- + q : ndarray + Charges to evaluate the requirement curve at + """ + req = self.requirement(q) + self.ax.plot(q, req, '--', color='black', label="Requirement") + + def plot_poisson(self, q): + """ + Plot the Poisson limit curve onto the figure. + + Parameters + ---------- + q : ndarray + Charges to evaluate the limit at + """ + poisson = self.poisson(q) + self.ax.plot(q, poisson, '--', color='grey', label="Poisson") + + +class ChargeResolutionWRRPlotter(ChargeResolutionPlotter): + def __init__(self, config=None, tool=None, **kwargs): + """ + Plots the charge resolution similarly to ChargeResolutionPlotter, with + the values divided by the CTA requirement. + + Parameters + ---------- + config : traitlets.loader.Config + Configuration specified by config file or cmdline arguments. + Used to set traitlet values. + Set to None if no configuration to pass. + tool : ctapipe.core.Tool + Tool executable that is calling this component. + Passes the correct logger to the component. + Set to None if no Tool to pass. + kwargs + """ + super().__init__(config=config, tool=tool, **kwargs) + self.ax.set_xlabel("True Charge (p.e.)") + self.ax.set_ylabel(r"$\frac{{\sigma_Q}}{{Q}}$ / Requirement") + + def _plot(self, x, y, **kwargs): + y = y / self.requirement(x) + if 'yerr' in kwargs: + kwargs['yerr'] /= self.requirement(x) + super()._plot(x, y, **kwargs) + + def plot_requirement(self, q): + req = self.requirement(q) + req /= self.requirement(q) + self.ax.plot(q, req, '--', color='black', label="Requirement") + + def plot_poisson(self, q): + poisson = self.poisson(q) + poisson /= self.requirement(q) + self.ax.plot(q, poisson, '--', color='grey', label="Poisson") + + def _finish(self): + self.ax.set_xscale('log') + self.ax.get_xaxis().set_major_formatter( + FuncFormatter(lambda x, _: '{:g}'.format(x))) diff --git a/ctapipe/plotting/tests/test_charge_resolution.py b/ctapipe/plotting/tests/test_charge_resolution.py new file mode 100644 index 00000000000..4d6cc063105 --- /dev/null +++ b/ctapipe/plotting/tests/test_charge_resolution.py @@ -0,0 +1,124 @@ +from ctapipe.analysis.camera.charge_resolution import \ + ChargeResolutionCalculator +import os +import numpy as np +import pandas as pd +from ..charge_resolution import ( + sum_errors, bin_dataframe, + ChargeResolutionPlotter, ChargeResolutionWRRPlotter +) +from numpy.testing import assert_almost_equal +import pytest + + +def create_temp_cr_file(directory): + chargeres = ChargeResolutionCalculator() + true_charge = np.arange(100) + measured_charge = np.arange(100) + chargeres.add(0, true_charge, measured_charge) + chargeres.add(0, true_charge, measured_charge) + df_p, df_c = chargeres.finish() + + output_path = os.path.join(str(directory), "cr.h5") + with pd.HDFStore(output_path, 'w') as store: + store['charge_resolution_pixel'] = df_p + store['charge_resolution_camera'] = df_c + return output_path + + +def test_sum_errors(): + errors = np.array([2, 5, 6, 7]) + assert_almost_equal(sum_errors(errors), 5.339, 3) + + +def test_bin_dataframe(): + chargeres = ChargeResolutionCalculator() + true_charge = np.arange(100) + measured_charge = np.arange(100) + chargeres.add(0, true_charge, measured_charge) + chargeres.add(0, true_charge, measured_charge) + df_p, df_c = chargeres.finish() + + df = bin_dataframe(df_p, 20) + assert 'bin' in df.columns + assert np.unique(df['bin']).size <= 20 + + +def test_file_reading(tmpdir): + path = create_temp_cr_file(tmpdir) + output_path = os.path.join(str(tmpdir), "cr.pdf") + plotter = ChargeResolutionPlotter(output_path=output_path) + plotter._set_file(path) + assert plotter._df_pixel is not None + assert plotter._df_camera is not None + + +def test_incorrect_file(tmpdir): + path = os.path.join(str(tmpdir), "cr_incorrect.h5") + output_path = os.path.join(str(tmpdir), "cr_incorrect.pdf") + with pd.HDFStore(path, 'w') as store: + store['test'] = pd.DataFrame(dict(a=[3])) + + plotter = ChargeResolutionPlotter(output_path=output_path) + with pytest.raises(KeyError): + plotter._set_file(path) + + +def test_missing_file(tmpdir): + path = os.path.join(str(tmpdir), "cr_missing.h5") + output_path = os.path.join(str(tmpdir), "cr_missing.pdf") + + assert not os.path.exists(path) + + plotter = ChargeResolutionPlotter(output_path=output_path) + with pytest.raises(OSError): + plotter._set_file(path) + + +def test_plotting(tmpdir): + path = create_temp_cr_file(tmpdir) + output_path = os.path.join(str(tmpdir), "cr.pdf") + plotter = ChargeResolutionPlotter(output_path=output_path) + plotter.plot_average(path, "average") + plotter.plot_camera(path, "average") + plotter.plot_pixel(path, 0, "average") + q = np.arange(1, 1000) + plotter.plot_requirement(q) + plotter.plot_poisson(q) + + plotter.save() + assert os.path.exists(output_path) + + +def test_limit_curves(): + value = ChargeResolutionPlotter.limit_curves(1, 2, 3, 4, 5, 6) + assert_almost_equal(value, 9.798, 3) + + +def test_requirement(): + value = ChargeResolutionPlotter.requirement(np.arange(1, 10))[0] + value2 = ChargeResolutionPlotter.requirement(np.arange(1, 10))[-1] + assert_almost_equal(value, 2.020, 3) + assert_almost_equal(value2, 0.450, 3) + + +def test_poisson(): + value = ChargeResolutionPlotter.poisson(np.arange(1, 10))[0] + value2 = ChargeResolutionPlotter.poisson(np.arange(1, 10))[-1] + assert_almost_equal(value, 1, 3) + assert_almost_equal(value2, 0.333, 3) + + +def test_plotting_wrr(tmpdir): + path = create_temp_cr_file(tmpdir) + output_path = os.path.join(str(tmpdir), "cr_wrr.pdf") + plotter = ChargeResolutionWRRPlotter(output_path=output_path) + plotter.plot_average(path, "average") + plotter.plot_camera(path, "average") + plotter.plot_pixel(path, 0, "average") + q = np.arange(1, 1000) + plotter.plot_requirement(q) + plotter.plot_poisson(q) + + plotter.save() + assert os.path.exists(output_path) diff --git a/ctapipe/tools/extract_charge_resolution.py b/ctapipe/tools/extract_charge_resolution.py index a578cdb2e74..f754756c116 100644 --- a/ctapipe/tools/extract_charge_resolution.py +++ b/ctapipe/tools/extract_charge_resolution.py @@ -1,54 +1,55 @@ """ -Extract data necessary to calcualte charge resolution from raw data files. +Calculate the Charge Resolution from a sim_telarray simulation and store +within a HDF5 file. """ import os - import numpy as np +import pandas as pd from tqdm import tqdm -from traitlets import Dict, List, Int, Unicode - -from ctapipe.analysis.camera.chargeresolution import ChargeResolutionCalculator +from traitlets import Dict, List, Unicode +from ctapipe.analysis.camera.charge_resolution import \ + ChargeResolutionCalculator from ctapipe.calib.camera.dl0 import CameraDL0Reducer from ctapipe.calib.camera.dl1 import CameraDL1Calibrator from ctapipe.calib.camera.r1 import HESSIOR1Calibrator -from ctapipe.core import Tool +from ctapipe.core import Tool, Provenance from ctapipe.image.charge_extractors import ChargeExtractorFactory from ctapipe.io.simteleventsource import SimTelEventSource class ChargeResolutionGenerator(Tool): name = "ChargeResolutionGenerator" - description = "Generate the a pickle file of ChargeResolutionFile for " \ - "a MC file." - - telescopes = List(Int, None, allow_none=True, - help='Telescopes to include from the event file. ' - 'Default = All telescopes').tag(config=True) - output_name = Unicode('charge_resolution', - help='Name of the output charge resolution hdf5 ' - 'file').tag(config=True) - - aliases = Dict(dict(f='SimTelEventSource.input_url', - max_events='SimTelEventSource.max_events', - extractor='ChargeExtractorFactory.product', - window_width='ChargeExtractorFactory.window_width', - t0='ChargeExtractorFactory.t0', - window_shift='ChargeExtractorFactory.window_shift', - sig_amp_cut_HG='ChargeExtractorFactory.sig_amp_cut_HG', - sig_amp_cut_LG='ChargeExtractorFactory.sig_amp_cut_LG', - lwt='ChargeExtractorFactory.lwt', - clip_amplitude='CameraDL1Calibrator.clip_amplitude', - radius='CameraDL1Calibrator.radius', - max_pe='ChargeResolutionCalculator.max_pe', - T='ChargeResolutionGenerator.telescopes', - O='ChargeResolutionGenerator.output_name', - )) - classes = List([SimTelEventSource, - ChargeExtractorFactory, - CameraDL1Calibrator, - ChargeResolutionCalculator - ]) + description = ("Calculate the Charge Resolution from a sim_telarray " + "simulation and store within a HDF5 file.") + + output_path = Unicode( + 'charge_resolution.h5', + help='Path to store the output HDF5 file' + ).tag(config=True) + + aliases = Dict(dict( + f='SimTelEventSource.input_url', + max_events='SimTelEventSource.max_events', + T='SimTelEventSource.allowed_tels', + extractor='ChargeExtractorFactory.product', + window_width='ChargeExtractorFactory.window_width', + t0='ChargeExtractorFactory.t0', + window_shift='ChargeExtractorFactory.window_shift', + sig_amp_cut_HG='ChargeExtractorFactory.sig_amp_cut_HG', + sig_amp_cut_LG='ChargeExtractorFactory.sig_amp_cut_LG', + lwt='ChargeExtractorFactory.lwt', + clip_amplitude='CameraDL1Calibrator.clip_amplitude', + radius='CameraDL1Calibrator.radius', + max_pe='ChargeResolutionCalculator.max_pe', + o='ChargeResolutionGenerator.output_path', + )) + classes = List([ + SimTelEventSource, + ChargeExtractorFactory, + CameraDL1Calibrator, + ChargeResolutionCalculator + ]) def __init__(self, **kwargs): super().__init__(**kwargs) @@ -72,50 +73,48 @@ def setup(self): self.dl1 = CameraDL1Calibrator(extractor=extractor, **kwargs) - self.calculator = ChargeResolutionCalculator(**kwargs) + self.calculator = ChargeResolutionCalculator() def start(self): - desc = "Filling Charge Resolution" + desc = "Extracting Charge Resolution" for event in tqdm(self.eventsource, desc=desc): - tels = list(event.dl0.tels_with_data) + self.r1.calibrate(event) + self.dl0.reduce(event) + self.dl1.calibrate(event) # Check events have true charge included if event.count == 0: try: - if np.all(event.mc.tel[ - tels[0]].photo_electron_image == 0): + pe = list(event.mc.tel.values())[0].photo_electron_image + if np.all(pe == 0): raise KeyError except KeyError: - self.log.exception('Source does not contain ' - 'true charge!') + self.log.exception( + 'Source does not contain true charge!' + ) raise - self.r1.calibrate(event) - self.dl0.reduce(event) - self.dl1.calibrate(event) - - if self.telescopes: - tels = [] - for tel in self.telescopes: - if tel in event.dl0.tels_with_data: - tels.append(tel) - - for telid in tels: - true_charge = event.mc.tel[telid].photo_electron_image - measured_charge = event.dl1.tel[telid].image[0] - self.calculator.add_charges(true_charge, measured_charge) + for mc, dl1 in zip(event.mc.tel.values(), event.dl1.tel.values()): + true_charge = mc.photo_electron_image + measured_charge = dl1.image[0] + pixels = np.arange(measured_charge.size) + self.calculator.add(pixels, true_charge, measured_charge) def finish(self): - input_url = self.eventsource.input_url - input_directory = os.path.dirname(input_url) - input_name = os.path.splitext(os.path.basename(input_url))[0] - output_directory = os.path.join(input_directory, input_name) + df_p, df_c = self.calculator.finish() + + output_directory = os.path.dirname(self.output_path) if not os.path.exists(output_directory): self.log.info("Creating directory: {}".format(output_directory)) os.makedirs(output_directory) - name = "{}.h5".format(self.output_name) - ouput_path = os.path.join(output_directory, name) - self.calculator.save(ouput_path) + + with pd.HDFStore(self.output_path, 'w') as store: + store['charge_resolution_pixel'] = df_p + store['charge_resolution_camera'] = df_c + + self.log.info("Created charge resolution file: {}" + .format(self.output_path)) + Provenance().add_output_file(self.output_path) def main(): diff --git a/ctapipe/tools/plot_charge_resolution.py b/ctapipe/tools/plot_charge_resolution.py index 9f749d70548..da2141cc6ff 100644 --- a/ctapipe/tools/plot_charge_resolution.py +++ b/ctapipe/tools/plot_charge_resolution.py @@ -1,172 +1,30 @@ """ -Plot charge resolution data generated by ctapipe-chargeres-extract. +Plot charge resolutions generated by ChargeResolutionCalculator. """ - -from os.path import dirname, basename, expanduser, exists, splitext -from os import makedirs -from matplotlib import pyplot as plt -from math import log10 -import warnings import numpy as np -from traitlets import Dict, List, Unicode, Int, Bool -from ctapipe.core import Tool, Component -from ctapipe.analysis.camera.chargeresolution import ChargeResolutionCalculator - - -class ChargeResolutionPlotter(Component): - output_path = Unicode(None, allow_none=True, - help='Output path to save the ' - 'plot.').tag(config=True) - max_pe = Int(2000, help='Maximum pe to plot').tag(config=True) - linear_x = Bool(False, help='Plot the x values on a linear axis, ' - 'instead of log').tag(config=True) - linear_y = Bool(False, help='Plot the y values on a linear axis, ' - 'instead of log').tag(config=True) - - def __init__(self, config=None, tool=None, **kwargs): - """ - Calculator of charge resolution. - - Parameters - ---------- - config : traitlets.loader.Config - Configuration specified by config file or cmdline arguments. - Used to set traitlet values. - Set to None if no configuration to pass. - tool : ctapipe.core.Tool - Tool executable that is calling this component. - Passes the correct logger to the component. - Set to None if no Tool to pass. - kwargs - """ - super().__init__(config=config, parent=tool, **kwargs) - - try: - if self.output_path is None: - raise ValueError - except ValueError: - self.log.exception('Please specify an output path') - raise - - self.fig = plt.figure(figsize=(20, 8)) - self.ax_l = self.fig.add_subplot(121) - self.ax_r = self.fig.add_subplot(122) - - self.fig.subplots_adjust(left=0.05, right=0.95, wspace=0.6) - - self.legend_handles = [] - self.legend_labels = [] - - def plot_chargeres(self, name, x, res, error): - valid = (x <= self.max_pe) - x_v = x[valid] - res_v = res[valid] - error_v = error[valid] - - eb_l, _, _ = self.ax_l.errorbar(x_v, res_v, yerr=error_v, - marker='x', linestyle="None") - self.legend_handles.append(eb_l) - self.legend_labels.append(splitext(name)[0]) - - def plot_scaled_chargeres(self, x, res, error): - valid = (x <= self.max_pe) - x_v = x[valid] - res_v = res[valid] - error_v = error[valid] - - self.ax_r.errorbar(x_v, res_v, yerr=error_v, - marker='x', linestyle="None") - - def plot_limit_curves(self): - x = np.logspace(log10(0.9), log10(self.max_pe * 1.1), 100) - requirement = ChargeResolutionCalculator.requirement(x) - goal = ChargeResolutionCalculator.goal(x) - poisson = ChargeResolutionCalculator.poisson(x) - - r_p, = self.ax_l.plot(x, requirement, 'r', ls='--') - g_p, = self.ax_l.plot(x, goal, 'g', ls='--') - p_p, = self.ax_l.plot(x, poisson, c='0.75', ls='--') - self.ax_r.plot(x, requirement / goal, 'r') - self.ax_r.plot(x, np.ones_like(x), 'g') - self.ax_r.plot(x, poisson / goal, c='0.75', ls='--') - - self.legend_handles.append(r_p) - self.legend_labels.append("Requirement") - self.legend_handles.append(g_p) - self.legend_labels.append("Goal") - self.legend_handles.append(p_p) - self.legend_labels.append("Poisson") - - def save(self): - if not self.linear_x: - self.ax_l.set_xscale('log') - self.ax_r.set_xscale('log') - if not self.linear_y: - self.ax_l.set_yscale('log') - - self.ax_l.set_xlabel(r'True Charge $Q_T$ (p.e.)') - self.ax_l.set_ylabel('Charge Resolution') - self.ax_r.set_xlabel(r'True Charge $Q_T$ (p.e.)') - self.ax_r.set_ylabel('Charge Resolution/Goal') - - self.ax_l.legend(self.legend_handles, self.legend_labels, - bbox_to_anchor=(1.02, 1.), loc=2, - borderaxespad=0., fontsize=10) - - self.fig.suptitle(splitext(basename(self.output_path))[0]) - - self.ax_l.set_xlim(0.9, self.max_pe * 1.1) - self.ax_r.set_xlim(0.9, self.max_pe * 1.1) - if self.max_pe > 2000: - self.ax_r.set_xlim(0.9, 2000 * 1.1) - # if args.maxpeplot is not None: - # ax_l.set_xlim(0.9, args.maxpeplot) - # ax_r.set_xlim(0.9, args.maxpeplot) - - warnings.filterwarnings("ignore", module="matplotlib") - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - if self.output_path is None: - plt.show() - else: - self.output_path = expanduser(self.output_path) - output_dir = dirname(self.output_path) - if not exists(output_dir): - self.log.info("[output] Creating directory: " - "{}".format(output_dir)) - makedirs(output_dir) - self.log.info("[output] {}".format(self.output_path)) - self.fig.savefig(self.output_path, bbox_inches='tight') +from traitlets import Dict, List, Unicode +from ctapipe.core import Tool +from ctapipe.plotting.charge_resolution import ChargeResolutionPlotter class ChargeResolutionViewer(Tool): name = "ChargeResolutionViewer" - description = "Plot the charge resolution from " \ - "ChargeResolutionCalculator objects restored via " \ - "pickled dictionaries." - - input_files = List(Unicode, None, - help='Input hdf5 files that are produced from ' - 'ChargeResolutionCalculator.save().' - '').tag(config=True) - - aliases = Dict(dict(f='ChargeResolutionViewer.input_files', - B='ChargeResolutionCalculator.binning', - max_pe='ChargeResolutionPlotter.max_pe', - O='ChargeResolutionPlotter.output_path', - )) - flags = Dict(dict(L=({'ChargeResolutionCalculator': {'log_bins': False}}, - 'Bin the x axis linearly instead of logarithmic.'), - linx=({'ChargeResolutionPlotter': {'linear_x': True}}, - 'Plot the x values on a linear axis, ' - 'instead of log.'), - liny=({'ChargeResolutionPlotter': {'linear_y': True}}, - 'Plot the x values on a linear axis, ' - 'instead of log.') - )) - classes = List([ChargeResolutionCalculator, - ChargeResolutionPlotter - ]) + description = ("Plot charge resolutions generated by " + "ChargeResolutionCalculator.") + + input_files = List( + Unicode, None, + help='Input HDF5 files produced by ChargeResolutionCalculator' + ).tag(config=True) + + aliases = Dict(dict( + f='ChargeResolutionViewer.input_files', + B='ChargeResolutionPlotter.n_bins', + o='ChargeResolutionPlotter.output_path', + )) + classes = List([ + ChargeResolutionPlotter, + ]) def __init__(self, **kwargs): super().__init__(**kwargs) @@ -177,21 +35,16 @@ def setup(self): self.log_format = "%(levelname)s: %(message)s [%(name)s.%(funcName)s]" kwargs = dict(config=self.config, tool=self) - self.calculator = ChargeResolutionCalculator(**kwargs) self.plotter = ChargeResolutionPlotter(**kwargs) def start(self): - self.plotter.plot_limit_curves() for fp in self.input_files: - self.calculator.load(fp) - x, res, res_error, scaled_res, scaled_res_error = \ - self.calculator.get_charge_resolution() - - name = basename(fp) - self.plotter.plot_chargeres(name, x, res, res_error) - self.plotter.plot_scaled_chargeres(x, scaled_res, scaled_res_error) + self.plotter.plot_camera(fp) def finish(self): + q = np.arange(1, 1000) + self.plotter.plot_poisson(q) + self.plotter.plot_requirement(q) self.plotter.save() diff --git a/ctapipe/tools/plot_charge_resolution_variation_hist.py b/ctapipe/tools/plot_charge_resolution_variation_hist.py deleted file mode 100644 index 2a8b5a3bdca..00000000000 --- a/ctapipe/tools/plot_charge_resolution_variation_hist.py +++ /dev/null @@ -1,146 +0,0 @@ -from os.path import dirname, exists, splitext, basename -from os import makedirs -from math import ceil, floor -from matplotlib import pyplot as plt -from math import log10 -import warnings -import numpy as np -from matplotlib.colors import LogNorm -from traitlets import Dict, List, Unicode -from ctapipe.core import Tool, Component -from ctapipe.analysis.camera.chargeresolution import ChargeResolutionCalculator - - -class ChargeResolutionVariationPlotter(Component): - output_path = Unicode(None, allow_none=True, - help='Output path to save the ' - 'plot.').tag(config=True) - - def __init__(self, config=None, tool=None, **kwargs): - """ - Calculator of charge resolution. - - Parameters - ---------- - config : traitlets.loader.Config - Configuration specified by config file or cmdline arguments. - Used to set traitlet values. - Set to None if no configuration to pass. - tool : ctapipe.core.Tool - Tool executable that is calling this component. - Passes the correct logger to the component. - Set to None if no Tool to pass. - kwargs - """ - super().__init__(config=config, parent=tool, **kwargs) - - try: - if self.output_path is None: - raise ValueError - except ValueError: - self.log.exception('Please specify an output path') - raise - - self.fig = plt.figure(figsize=(20, 8)) - self.ax_l = self.fig.add_subplot(121) - self.ax_r = self.fig.add_subplot(122) - - self.fig.subplots_adjust(left=0.05, right=0.95, wspace=0.6) - - self.legend_handles = [] - self.legend_labels = [] - - def plot_hist(self, hist, xedges, yedges): - - hist[hist == 0.0] = np.nan - - fig = plt.figure(figsize=(10, 7)) - ax = fig.add_subplot(111) - ax.set_title(splitext(basename(self.output_path))[0]) - x, y = np.meshgrid(xedges, yedges) - x = np.power(10, x) - y = np.power(10, y) - hist_mask = np.ma.masked_where(np.isnan(hist), hist) - im = ax.pcolormesh(x, y, hist_mask, norm=LogNorm()) - cb = plt.colorbar(im) - ax.set_aspect('equal') - ax.grid() - ax.set_xscale('log') - ax.set_yscale('log') - ax.set_xlabel(r'True Charge $Q_T$ (p.e.)') - ax.set_ylabel(r'Measured Charge $Q_M$ (p.e.)') - cb.ax.set_ylabel("Count") - - line = np.linspace(*ax.get_xlim(), 100) - ax.plot(line, line, c='0.75', ls='--') - - # Add minor ticks - lmin = floor(log10(hist_mask.min())) - lmax = ceil(log10(hist_mask.max())) - logticks = np.tile(np.arange(lmin, 10), lmax) * ( - np.power(10, np.arange(lmax * 10) // 10)) - logticks = im.norm(logticks[(logticks != 0) & - (logticks >= hist_mask.min()) & - (logticks <= hist_mask.max())]) - cb.ax.yaxis.set_ticks(logticks, minor=True) - cb.ax.yaxis.set_tick_params(which='minor', length=5) - cb.ax.tick_params(length=10) - - output_dir = dirname(self.output_path) - if not exists(output_dir): - self.log.info("[output] Creating directory: {}".format(output_dir)) - makedirs(output_dir) - self.log.info("[output] {}".format(self.output_path)) - warnings.filterwarnings("ignore", module="matplotlib") - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - fig.savefig(self.output_path, bbox_inches='tight') - - -class ChargeResolutionVariationViewer(Tool): - name = "ChargeResolutionVariationViewer" - description = "Plot the charge resolution from " \ - "ChargeResolutionCalculator objects restored via " \ - "pickled dictionaries." - - input_path = Unicode(None, allow_none=True, - help='Path to the hdf5 file produced from' - 'ChargeResolutionCalculator.save()' - '').tag(config=True) - - aliases = Dict(dict(f='ChargeResolutionVariationViewer.input_path', - O='ChargeResolutionVariationPlotter.output_path', - )) - classes = List([ChargeResolutionCalculator, - ChargeResolutionVariationPlotter - ]) - - def __init__(self, **kwargs): - super().__init__(**kwargs) - self.calculator = None - self.plotter = None - - def setup(self): - self.log_format = "%(levelname)s: %(message)s [%(name)s.%(funcName)s]" - kwargs = dict(config=self.config, tool=self) - - self.calculator = ChargeResolutionCalculator(**kwargs) - self.plotter = ChargeResolutionVariationPlotter(**kwargs) - - def start(self): - self.calculator.load(self.input_path) - - def finish(self): - hist = self.calculator.variation_hist - xedges = self.calculator.variation_xedges - yedges = self.calculator.variation_yedges - self.plotter.plot_hist(hist, xedges, yedges) - - -def main(): - exe = ChargeResolutionVariationViewer() - exe.run() - - -if __name__ == '__main__': - main() diff --git a/ctapipe/tools/tests/test_tools.py b/ctapipe/tools/tests/test_tools.py index 1ab7dcdf6db..fe0f252799a 100644 --- a/ctapipe/tools/tests/test_tools.py +++ b/ctapipe/tools/tests/test_tools.py @@ -3,7 +3,11 @@ from ctapipe.tools.dump_instrument import DumpInstrumentTool from ctapipe.tools.info import info from ctapipe.tools.bokeh.file_viewer import BokehFileViewer +from ctapipe.tools.extract_charge_resolution import ChargeResolutionGenerator +from ctapipe.tools.plot_charge_resolution import ChargeResolutionViewer from ctapipe.utils import get_dataset_path +import os +import pytest def test_info(): @@ -49,3 +53,29 @@ def test_bokeh_file_viewer(): tool.run() assert tool.reader.input_url == get_dataset_path("gamma_test.simtel.gz") + + +def test_extract_charge_resolution(tmpdir): + output_path = os.path.join(str(tmpdir), "cr.h5") + tool = ChargeResolutionGenerator() + with pytest.raises(KeyError): + tool.run([ + '-f', get_dataset_path("gamma_test_large.simtel.gz"), + '-o', output_path, + ]) + # TODO: Test files do not contain true charge, cannot test tool fully + # assert os.path.exists(output_path) + + +def test_plot_charge_resolution(tmpdir): + from ctapipe.plotting.tests.test_charge_resolution import \ + create_temp_cr_file + path = create_temp_cr_file(tmpdir) + + output_path = os.path.join(str(tmpdir), "cr.pdf") + tool = ChargeResolutionViewer() + tool.run([ + '-f', [path], + '-o', output_path, + ]) + assert os.path.exists(output_path) diff --git a/setup.py b/setup.py index 9edb5606397..6524fcc8cea 100755 --- a/setup.py +++ b/setup.py @@ -33,8 +33,6 @@ 'ctapipe-dump-triggers = ctapipe.tools.dump_triggers:main', 'ctapipe-chargeres-extract = ctapipe.tools.extract_charge_resolution:main', 'ctapipe-chargeres-plot = ctapipe.tools.plot_charge_resolution:main', - 'ctapipe-chargeres-hist = ' - 'ctapipe.tools.plot_charge_resolution_variation_hist:main', 'ctapipe-dump-instrument=ctapipe.tools.dump_instrument:main', 'ctapipe-event-viewer = ctapipe.tools.bokeh.file_viewer:main' ] From bc12d6bbdcecfe8551124356a38d3a723f302558 Mon Sep 17 00:00:00 2001 From: Dominik Neise Date: Mon, 21 Jan 2019 11:37:49 +0100 Subject: [PATCH 4/7] remove optional deps in order to fix tests (#925) * remove optional deps in order to fix tests * remove tests for py3.5 since it does actually test for 3.6 only * add py3.7 and osx to test matrix * xvfb does not work on osx * allow this to take 30min * fix typo its travis_wait .. not what I wrote before * remove tests for code, which we do not even want anymore * modify the default to SimTelEventSource * re-add line accidentally removed before * install pyzmq>=17 needed for notebooks on OSX * Add env files for each python version * Pip install pyhessio for now * Fix pyhessio url * Fix pyhessio version * remove osx tests for the time being * Revert "remove tests for code, which we do not even want anymore" This reverts commit 766f5e213db42811045684cec38e797341e64f2f. --- .travis.yml | 58 +++++++++++++-------------------- environment.yml | 1 + examples/simple_event_writer.py | 2 +- py3.6_env.yaml | 46 ++++++++++++++++++++++++++ py3.7_env.yaml | 46 ++++++++++++++++++++++++++ 5 files changed, 116 insertions(+), 37 deletions(-) create mode 100644 py3.6_env.yaml create mode 100644 py3.7_env.yaml diff --git a/.travis.yml b/.travis.yml index 82f87bce8a8..a99a4f62db2 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,21 +1,24 @@ -language: python -sudo: required +language: generic -python: - # not used directly, but this sets TRAVIS_PYTHON_VERSION so we can use it - # in anaconda as well (we don't support python less than 3.5) - - 3.5 - - 3.6 +matrix: + include: + - os: linux + python: 3.6 + env: + - PYTHON_VERSION=3.6 + - os: linux + python: 3.7 + env: + - PYTHON_VERSION=3.7 + # - os: osx + # python: 3.6 + # env: + # - PYTHON_VERSION=3.6 + # - os: osx + # python: 3.7 + # env: + # - PYTHON_VERSION=3.7 -os: - - linux - # - osx # currently osx python projects are not supported in Travis - -env: - global: - - TARGET_SOFTWARE=$HOME/Software/target_software - - LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$HOME/Software/TargetDriver/install/lib:$HOME/Software/TargetIO/install/lib:$HOME/Software/TargetCalib/install/lib - - TARGETCALIBPATH=$HOME/Software/TargetCalib/install before_install: @@ -23,7 +26,7 @@ before_install: # against future changes - export PYTHONIOENCODING=UTF8 - - export MPLBACKEND=agg + - export MPLBACKEND=Agg # Install miniconda following instructions at # http://conda.pydata.org/docs/travis.html @@ -40,32 +43,15 @@ before_install: - conda update -q conda # get latest conda version # Useful for debugging any issues with conda - conda info -a - - # Make sure that interactive matplotlib backends work - - export DISPLAY=:99.0 - - sh -e /etc/init.d/xvfb start - git fetch --tags install: - - conda create --name cta-dev python=$TRAVIS_PYTHON_VERSION - - conda env update -n cta-dev --file environment.yml + - conda create --name cta-dev python=$PYTHON_VERSION + - travis_wait 20 conda env update -n cta-dev --file py${PYTHON_VERSION}_env.yaml - source activate cta-dev - ulimit -s 16000 # increase stack size limit, for libhessio - pip install travis-sphinx - pip install codecov - # ----- SST1M: - - pip install https://github.com/cta-sst-1m/protozfitsreader/archive/v1.4.2.tar.gz - # ----- end of SST1M - # ----- target_software - - conda install -c conda-forge cfitsio - - conda install swig - - mkdir -p $HOME/Software - - cd $HOME/Software - - git clone https://github.com/watsonjj/target_software.git - - cd $TARGET_SOFTWARE - - ./install.sh - - cd $TRAVIS_BUILD_DIR - # ----- end of target_software - python setup.py develop script: diff --git a/environment.yml b/environment.yml index 087142e26b8..10c54ac0c87 100644 --- a/environment.yml +++ b/environment.yml @@ -5,6 +5,7 @@ channels: - conda-forge dependencies: - zeromq + - pyzmq>=17 # needed for correct function of notebooks on OSX - astropy - bokeh - cython diff --git a/examples/simple_event_writer.py b/examples/simple_event_writer.py index 8cd2f2ed79b..09fd4f927ba 100755 --- a/examples/simple_event_writer.py +++ b/examples/simple_event_writer.py @@ -38,7 +38,7 @@ def setup(self): self.log.info('Configure EventSourceFactory...') self.event_source = EventSourceFactory.produce( - config=self.config, tool=self, product='HESSIOEventSource' + config=self.config, tool=self, product='SimTelEventSource' ) self.event_source.allowed_tels = self.config['Analysis']['allowed_tels'] diff --git a/py3.6_env.yaml b/py3.6_env.yaml new file mode 100644 index 00000000000..68d905a2e9d --- /dev/null +++ b/py3.6_env.yaml @@ -0,0 +1,46 @@ +name: cta-dev +channels: + - default + - cta-observatory + - conda-forge +dependencies: + - python=3.6 # nail the python version, so conda does not try upgrading / dowgrading + - zeromq + - pyzmq>=17 # needed for correct function of notebooks on OSX + - astropy + - bokeh + - cython + - gammapy + - graphviz + - h5py + - iminuit + - ipython + - ipywidgets + - jupyter + - matplotlib + - numba + - numpy>=1.15.4 + - numpydoc + - pandas + - pytest + - pytest-cov + - psutil + - scikit-learn + - scipy + - setuptools + - sphinx=1.7 + - sphinx_rtd_theme + - sphinx-automodapi + - traitlets + - wheel + - xz + - pyyaml + - zlib + - iminuit + - pytables + - tqdm + - pip: + - pytest_runner + - eventio==0.11.0 + - https://github.com/cta-observatory/ctapipe-extra/archive/v0.2.16.tar.gz + - https://github.com/cta-observatory/pyhessio/archive/v2.1.1.tar.gz diff --git a/py3.7_env.yaml b/py3.7_env.yaml new file mode 100644 index 00000000000..21f806a11d3 --- /dev/null +++ b/py3.7_env.yaml @@ -0,0 +1,46 @@ +name: cta-dev +channels: + - default + - cta-observatory + - conda-forge +dependencies: + - python=3.7 + - zeromq + - pyzmq>=17 # needed for correct function of notebooks on OSX + - astropy + - bokeh + - cython + - gammapy + - graphviz + - h5py + - iminuit + - ipython + - ipywidgets + - jupyter + - matplotlib + - numba + - numpy>=1.15.4 + - numpydoc + - pandas + - pytest + - pytest-cov + - psutil + - scikit-learn + - scipy + - setuptools + - sphinx=1.7 + - sphinx_rtd_theme + - sphinx-automodapi + - traitlets + - wheel + - xz + - pyyaml + - zlib + - iminuit + - pytables + - tqdm + - pip: + - pytest_runner + - eventio==0.11.0 + - https://github.com/cta-observatory/ctapipe-extra/archive/v0.2.16.tar.gz + - https://github.com/cta-observatory/pyhessio/archive/v2.1.1.tar.gz From c2efb649a337e1898ae072614cb1c1af49ecf9fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Wed, 23 Jan 2019 09:40:31 +0100 Subject: [PATCH 5/7] Prefix containers (#833) * Give containers a prefix attribute * Allow setting a prefix for container columns in hdfwriter * Add support for instance level container prefixes --- ctapipe/core/container.py | 81 ++++++++------ ctapipe/core/tests/test_container.py | 35 ++++++ ctapipe/io/containers.py | 4 + ctapipe/io/hdf5tableio.py | 41 +++++-- ctapipe/io/tableio.py | 4 +- ctapipe/io/tests/test_hdf5.py | 151 +++++++++++++++----------- ctapipe/plotting/tests/test_camera.py | 28 +++-- 7 files changed, 219 insertions(+), 125 deletions(-) diff --git a/ctapipe/core/container.py b/ctapipe/core/container.py index 562aa71f81e..43606921e36 100644 --- a/ctapipe/core/container.py +++ b/ctapipe/core/container.py @@ -46,11 +46,11 @@ class ContainerMeta(type): and no new fields can be added to a container by accident. ''' def __new__(cls, name, bases, dct): - items = [ + field_names = [ k for k, v in dct.items() if isinstance(v, Field) ] - dct['__slots__'] = tuple(items + ['meta']) + dct['__slots__'] = tuple(field_names + ['meta', 'prefix']) dct['fields'] = {} # inherit fields from baseclasses @@ -59,10 +59,16 @@ def __new__(cls, name, bases, dct): for k, v in b.fields.items(): dct['fields'][k] = v - for k in items: + for k in field_names: dct['fields'][k] = dct.pop(k) - return type.__new__(cls, name, bases, dct) + new_cls = type.__new__(cls, name, bases, dct) + + # if prefix was not set as a class variable, build a default one + if 'container_prefix' not in dct: + new_cls.container_prefix = name.lower().replace('container', '') + + return new_cls class Container(metaclass=ContainerMeta): @@ -72,7 +78,7 @@ class Container(metaclass=ContainerMeta): The purpose of this class is to provide a flexible data structure that works a bit like a dict or blank Python class, but prevents the user from accessing members that have not been defined a - priori (more like a C struct), and also keeps metdata information + priori (more like a C struct), and also keeps metadata information such as a description, defaults, and units for each item in the container. @@ -114,8 +120,12 @@ class Container(metaclass=ContainerMeta): """ def __init__(self, **fields): - self.meta = {} + # __slots__ cannot be provided with defaults + # via class variables, so we use a `__prefix` class variable + # and a `_prefix` in `__slots__` together with a property. + self.prefix = self.container_prefix + for k, v in self.fields.items(): setattr(self, k, deepcopy(v.default)) @@ -128,9 +138,12 @@ def __getitem__(self, key): def __setitem__(self, key, value): return setattr(self, key, value) - def items(self): + def items(self, add_prefix=False): """Generator over (key, value) pairs for the items""" - return ((k, getattr(self, k)) for k in self.fields.keys()) + if not add_prefix or self.prefix == '': + return ((k, getattr(self, k)) for k in self.fields.keys()) + + return ((self.prefix + '_' + k, getattr(self, k)) for k in self.fields.keys()) def keys(self): """Get the keys of the container""" @@ -140,7 +153,7 @@ def values(self): """Get the keys of the container""" return (getattr(self, k) for k in self.fields.keys()) - def as_dict(self, recursive=False, flatten=False): + def as_dict(self, recursive=False, flatten=False, add_prefix=False): """ convert the `Container` into a dictionary @@ -153,32 +166,27 @@ def as_dict(self, recursive=False, flatten=False): by appending the sub-Container name. """ if not recursive: - return dict(self.items()) + return dict(self.items(add_prefix=add_prefix)) else: d = dict() - for key, val in self.items(): - if key.startswith("_"): - continue + for key, val in self.items(add_prefix=add_prefix): if isinstance(val, Container) or isinstance(val, Map): if flatten: - d.update({"{}_{}".format(key, k): v - for k, v in val.as_dict(recursive).items()}) + d.update({ + "{}_{}".format(key, k): v + for k, v in val.as_dict( + recursive, + add_prefix=add_prefix + ).items() + }) else: - d[key] = val.as_dict(recursive=recursive, - flatten=flatten) - continue - d[key] = val + d[key] = val.as_dict( + recursive=recursive, flatten=flatten, add_prefix=add_prefix + ) + else: + d[key] = val return d - @classmethod - def disable_attribute_check(cls): - """ - Globally turn off attribute checking for all Containers, - which provides a ~5-10x speed up for setting attributes. - This may be used e.g. after code is tested to speed up operation. - """ - cls.__setattr__ = object.__setattr__ - def reset(self, recursive=True): """ set all values back to their default values""" for name, value in self.fields.items(): @@ -219,7 +227,7 @@ class Map(defaultdict): by `tel_id` or algorithm name). """ - def as_dict(self, recursive=False, flatten=False): + def as_dict(self, recursive=False, flatten=False, add_prefix=False): if not recursive: return dict(self.items()) else: @@ -227,11 +235,18 @@ def as_dict(self, recursive=False, flatten=False): for key, val in self.items(): if isinstance(val, Container) or isinstance(val, Map): if flatten: - d.update({"{}_{}".format(key, k): v - for k, v in val.as_dict(recursive).items()}) + d.update({ + "{}_{}".format(key, k): v + for k, v in val.as_dict( + recursive, add_prefix=add_prefix + ).items() + }) else: - d[key] = val.as_dict(recursive=recursive, - flatten=flatten) + d[key] = val.as_dict( + recursive=recursive, + flatten=flatten, + add_prefix=add_prefix, + ) continue d[key] = val return d diff --git a/ctapipe/core/tests/test_container.py b/ctapipe/core/tests/test_container.py index 531f8e950b5..7f82f13dee0 100644 --- a/ctapipe/core/tests/test_container.py +++ b/ctapipe/core/tests/test_container.py @@ -3,6 +3,41 @@ from ctapipe.core import Container, Field, Map +def test_prefix(): + class AwesomeContainer(Container): + pass + + # make sure the default prefix is class name without container + assert AwesomeContainer.container_prefix == 'awesome' + assert AwesomeContainer().prefix == 'awesome' + + # make sure we can set the class level prefix at definition time + class ReallyAwesomeContainer(Container): + container_prefix = 'test' + + assert ReallyAwesomeContainer.container_prefix == 'test' + r = ReallyAwesomeContainer() + assert r.prefix == 'test' + + ReallyAwesomeContainer.container_prefix = 'test2' + # new instance should have the old prefix, old instance + # the one it was created with + assert ReallyAwesomeContainer().prefix == 'test2' + assert r.prefix == 'test' + + # Make sure we can set the class level prefix at runtime + ReallyAwesomeContainer.container_prefix = 'foo' + assert ReallyAwesomeContainer().prefix == 'foo' + + # make sure we can assign instance level prefixes + c1 = ReallyAwesomeContainer() + c2 = ReallyAwesomeContainer() + c2.prefix = 'c2' + + assert c1.prefix == 'foo' + assert c2.prefix == 'c2' + + def test_inheritance(): class ExampleContainer(Container): diff --git a/ctapipe/io/containers.py b/ctapipe/io/containers.py index ba6dcfe128e..0f04d975e2b 100644 --- a/ctapipe/io/containers.py +++ b/ctapipe/io/containers.py @@ -706,6 +706,8 @@ class MuonIntensityParameter(Container): class HillasParametersContainer(Container): + container_prefix = 'hillas' + intensity = Field(nan, 'total intensity (size)') x = Field(nan, 'centroid x coordinate') @@ -725,6 +727,8 @@ class LeakageContainer(Container): """ Leakage """ + container_prefix = '' + leakage1_pixel = Field( nan, 'Percentage of pixels after cleaning' diff --git a/ctapipe/io/hdf5tableio.py b/ctapipe/io/hdf5tableio.py index 3a5d4b17ac9..da36bdd8be8 100644 --- a/ctapipe/io/hdf5tableio.py +++ b/ctapipe/io/hdf5tableio.py @@ -60,6 +60,8 @@ class HDF5TableWriter(TableWriter): group_name: str name of group into which to put all of the tables generated by this Writer (it will be placed under "/" in the file) + add_prefix: bool + if True, add the container prefix before each column name mode : str ('w', 'a') 'w' if you want to overwrite the file 'a' if you want to append data to the file @@ -72,9 +74,17 @@ class HDF5TableWriter(TableWriter): """ - def __init__(self, filename, group_name, mode='w', root_uep='/', **kwargs): - - super().__init__() + def __init__( + self, + filename, + group_name, + add_prefix=False, + mode='w', + root_uep='/', + **kwargs + ): + + super().__init__(add_prefix=add_prefix) self._schemas = {} self._tables = {} @@ -129,18 +139,24 @@ class Schema(tables.IsDescription): # create pytables schema description for the given container for container in containers: - for col_name, value in container.items(): + for col_name, value in container.items(add_prefix=self.add_prefix): typename = "" shape = 1 if self._is_column_excluded(table_name, col_name): - self.log.debug("excluded column: %s/%s", - table_name, col_name) + self.log.debug( + "excluded column: %s/%s", table_name, col_name + ) continue if isinstance(value, Quantity): - req_unit = container.fields[col_name].unit + if self.add_prefix: + key = col_name.replace(container.prefix + '_', '') + else: + key = col_name + req_unit = container.fields[key].unit + if req_unit is not None: tr = partial(tr_convert_and_strip_unit, unit=req_unit) meta['{}_UNIT'.format(col_name)] = str(req_unit) @@ -207,14 +223,17 @@ def _append_row(self, table_name, containers): row = table.row for container in containers: - for colname in filter(lambda c: c in table.colnames, - container.keys()): + selected_fields = filter( + lambda kv: kv[0] in table.colnames, + container.items(add_prefix=self.add_prefix) + ) + for colname, value in selected_fields: + value = self._apply_col_transform( - table_name, colname, container[colname] + table_name, colname, value ) row[colname] = value - row.append() def write(self, table_name, containers): diff --git a/ctapipe/io/tableio.py b/ctapipe/io/tableio.py index 7d669b36cdd..89654a8d34c 100644 --- a/ctapipe/io/tableio.py +++ b/ctapipe/io/tableio.py @@ -18,17 +18,17 @@ class TableWriter(Component, metaclass=ABCMeta): - `ctapipe.io.HDF5TableWriter` """ - def __init__(self, parent=None, **kwargs): + def __init__(self, parent=None, add_prefix=False, **kwargs): super().__init__(parent, **kwargs) self._transforms = defaultdict(dict) self._exclusions = defaultdict(list) + self.add_prefix = add_prefix def __enter__(self): return self def __exit__(self, exc_type, exc_val, exc_tb): - self.close() def exclude(self, table_name, pattern): diff --git a/ctapipe/io/tests/test_hdf5.py b/ctapipe/io/tests/test_hdf5.py index 9b7f89886b5..b06295862de 100644 --- a/ctapipe/io/tests/test_hdf5.py +++ b/ctapipe/io/tests/test_hdf5.py @@ -3,10 +3,14 @@ import numpy as np import pytest import tables +import pandas as pd from astropy import units as u from ctapipe.core.container import Container, Field -from ctapipe.io.containers import R0CameraContainer, MCEventContainer +from ctapipe.io.containers import ( + R0CameraContainer, MCEventContainer, + HillasParametersContainer, LeakageContainer, +) from ctapipe.io.hdf5tableio import HDF5TableWriter, HDF5TableReader @@ -28,24 +32,49 @@ def test_write_container(temp_h5_file): r0tel.meta['test_attribute'] = 3.14159 r0tel.meta['date'] = "2020-10-10" - writer = HDF5TableWriter(str(temp_h5_file), group_name='R0', - filters=tables.Filters( - complevel=7)) - writer.exclude("tel_002", ".*samples") # test exclusion of columns - - for ii in range(100): - r0tel.waveform[:] = np.random.uniform(size=(50, 10)) - r0tel.image[:] = np.random.uniform(size=50) - r0tel.num_samples = 10 - mc.energy = 10**np.random.uniform(1, 2) * u.TeV - mc.core_x = np.random.uniform(-1, 1) * u.m - mc.core_y = np.random.uniform(-1, 1) * u.m - - writer.write("tel_001", r0tel) - writer.write("tel_002", r0tel) # write a second table too - writer.write("MC", mc) - - writer.close() + with HDF5TableWriter( + temp_h5_file, + group_name='R0', + filters=tables.Filters(complevel=7) + ) as writer: + writer.exclude("tel_002", ".*samples") # test exclusion of columns + + for ii in range(100): + r0tel.waveform[:] = np.random.uniform(size=(50, 10)) + r0tel.image[:] = np.random.uniform(size=50) + r0tel.num_samples = 10 + mc.energy = 10**np.random.uniform(1, 2) * u.TeV + mc.core_x = np.random.uniform(-1, 1) * u.m + mc.core_y = np.random.uniform(-1, 1) * u.m + + writer.write("tel_001", r0tel) + writer.write("tel_002", r0tel) # write a second table too + writer.write("MC", mc) + + +def test_prefix(tmp_path): + tmp_file = tmp_path / 'test_prefix.hdf5' + hillas_parameter_container = HillasParametersContainer( + x=1 * u.m, y=1 * u.m, length=1 * u.m, width=1 * u.m, + ) + + leakage_container = LeakageContainer( + leakage1_pixel=0.1, + leakage2_pixel=0.1, + leakage1_intensity=0.1, + leakage2_intensity=0.1, + ) + + with HDF5TableWriter( + tmp_file.name, + group_name='blabla', + add_prefix=True + ) as writer: + writer.write('events', [hillas_parameter_container, leakage_container]) + + df = pd.read_hdf(tmp_file.name) + assert 'hillas_x' in df.columns + assert 'leakage2_pixel' in df.columns def test_write_containers(temp_h5_file): @@ -59,16 +88,14 @@ class C2(Container): d = Field('d', None) with tempfile.NamedTemporaryFile() as f: - writer = HDF5TableWriter(f.name, 'test') - for i in range(20): - c1 = C1() - c2 = C2() - c1.a, c1.b, c2.c, c2.d = np.random.normal(size=4) - c1.b = np.random.normal() - - writer.write("tel_001", [c1, c2]) + with HDF5TableWriter(f.name, 'test') as writer: + for i in range(20): + c1 = C1() + c2 = C2() + c1.a, c1.b, c2.c, c2.d = np.random.normal(size=4) + c1.b = np.random.normal() - writer.close() + writer.write("tel_001", [c1, c2]) def test_read_container(temp_h5_file): @@ -76,41 +103,37 @@ def test_read_container(temp_h5_file): r0tel2 = R0CameraContainer() mc = MCEventContainer() - reader = HDF5TableReader(str(temp_h5_file)) - - # get the generators for each table - mctab = reader.read('/R0/MC', mc) - r0tab1 = reader.read('/R0/tel_001', r0tel1) - r0tab2 = reader.read('/R0/tel_002', r0tel2) + with HDF5TableReader(temp_h5_file) as reader: - # read all 3 tables in sync - for ii in range(3): + # get the generators for each table + mctab = reader.read('/R0/MC', mc) + r0tab1 = reader.read('/R0/tel_001', r0tel1) + r0tab2 = reader.read('/R0/tel_002', r0tel2) - m = next(mctab) - r0_1 = next(r0tab1) - r0_2 = next(r0tab2) + # read all 3 tables in sync + for ii in range(3): - print("MC:", m) - print("t0:", r0_1.image) - print("t1:", r0_2.image) - print("---------------------------") + m = next(mctab) + r0_1 = next(r0tab1) + r0_2 = next(r0tab2) - assert 'test_attribute' in r0_1.meta - assert r0_1.meta['date'] == "2020-10-10" + print("MC:", m) + print("t0:", r0_1.image) + print("t1:", r0_2.image) + print("---------------------------") - reader.close() + assert 'test_attribute' in r0_1.meta + assert r0_1.meta['date'] == "2020-10-10" def test_read_whole_table(temp_h5_file): mc = MCEventContainer() - reader = HDF5TableReader(str(temp_h5_file)) - - for cont in reader.read('/R0/MC', mc): - print(cont) + with HDF5TableReader(temp_h5_file) as reader: - reader.close() + for cont in reader.read('/R0/MC', mc): + print(cont) def test_with_context_writer(temp_h5_file): @@ -135,39 +158,41 @@ def test_writer_closes_file(temp_h5_file): with tempfile.NamedTemporaryFile() as f: with HDF5TableWriter(f.name, 'test') as h5_table: - assert h5_table._h5file.isopen == True + assert h5_table._h5file.isopen == 1 - assert h5_table._h5file.isopen == False + assert h5_table._h5file.isopen == 0 def test_reader_closes_file(temp_h5_file): - with HDF5TableReader(str(temp_h5_file)) as h5_table: + with HDF5TableReader(temp_h5_file) as h5_table: - assert h5_table._h5file.isopen == True + assert h5_table._h5file.isopen == 1 - assert h5_table._h5file.isopen == False + assert h5_table._h5file.isopen == 0 def test_with_context_reader(temp_h5_file): mc = MCEventContainer() - with HDF5TableReader(str(temp_h5_file)) as h5_table: + with HDF5TableReader(temp_h5_file) as h5_table: - assert h5_table._h5file.isopen == True + assert h5_table._h5file.isopen == 1 for cont in h5_table.read('/R0/MC', mc): print(cont) - assert h5_table._h5file.isopen == False + assert h5_table._h5file.isopen == 0 def test_closing_reader(temp_h5_file): - f = HDF5TableReader(str(temp_h5_file)) + f = HDF5TableReader(temp_h5_file) f.close() + assert f._h5file.isopen == 0 + def test_closing_writer(temp_h5_file): @@ -175,34 +200,32 @@ def test_closing_writer(temp_h5_file): h5_table = HDF5TableWriter(f.name, 'test') h5_table.close() + assert h5_table._h5file.isopen == 0 + def test_cannot_read_with_writer(temp_h5_file): with pytest.raises(IOError): with HDF5TableWriter(temp_h5_file, 'test', mode='r'): - pass def test_cannot_write_with_reader(temp_h5_file): with HDF5TableReader(temp_h5_file, mode='w') as h5: - assert h5._h5file.mode == 'r' def test_cannot_append_with_reader(temp_h5_file): with HDF5TableReader(temp_h5_file, mode='a') as h5: - assert h5._h5file.mode == 'r' def test_cannot_r_plus_with_reader(temp_h5_file): with HDF5TableReader(temp_h5_file, mode='r+') as h5: - assert h5._h5file.mode == 'r' diff --git a/ctapipe/plotting/tests/test_camera.py b/ctapipe/plotting/tests/test_camera.py index 2ca46ed43d8..10f4133e08e 100644 --- a/ctapipe/plotting/tests/test_camera.py +++ b/ctapipe/plotting/tests/test_camera.py @@ -6,26 +6,24 @@ def test_eventplotter(): dataset = get_dataset_path("gamma_test.simtel.gz") - with event_source(dataset, max_events=1) as source: event = next(iter(source)) - telid = list(event.r0.tels_with_data)[0] - - data = event.r0.tel[telid].waveform[0] - plotter = CameraPlotter(event) + telid = list(event.r0.tels_with_data)[0] - camera = plotter.draw_camera(telid, data[:, 0]) - assert camera is not None - np.testing.assert_array_equal(camera.image, data[:, 0]) + data = event.r0.tel[telid].waveform[0] + plotter = CameraPlotter(event) - plotter.draw_camera_pixel_ids(telid, [0, 1, 2]) + camera = plotter.draw_camera(telid, data[:, 0]) + assert camera is not None + np.testing.assert_array_equal(camera.image, data[:, 0]) - waveform = plotter.draw_waveform(data[0, :]) - assert waveform is not None - np.testing.assert_array_equal(waveform.get_ydata(), data[0, :]) + plotter.draw_camera_pixel_ids(telid, [0, 1, 2]) - line = plotter.draw_waveform_positionline(0) - assert line is not None - np.testing.assert_array_equal(line.get_xdata(), [0, 0]) + waveform = plotter.draw_waveform(data[0, :]) + assert waveform is not None + np.testing.assert_array_equal(waveform.get_ydata(), data[0, :]) + line = plotter.draw_waveform_positionline(0) + assert line is not None + np.testing.assert_array_equal(line.get_xdata(), [0, 0]) From 293fcd6b8a2ef7d821f8f390fada6e9cccb483c3 Mon Sep 17 00:00:00 2001 From: Dominik Neise Date: Wed, 23 Jan 2019 13:42:04 +0100 Subject: [PATCH 6/7] speed up is_compatible; makes event_source factory faster (#927) speed up is_compatible; makes event_source factory faster --- ctapipe/io/lsteventsource.py | 20 ++++++++++++-------- ctapipe/io/nectarcameventsource.py | 14 +++++++++----- ctapipe/io/sst1meventsource.py | 22 ++++++++++++++++++++++ 3 files changed, 43 insertions(+), 13 deletions(-) diff --git a/ctapipe/io/lsteventsource.py b/ctapipe/io/lsteventsource.py index a02e2f42b05..0246b5fa921 100644 --- a/ctapipe/io/lsteventsource.py +++ b/ctapipe/io/lsteventsource.py @@ -25,10 +25,10 @@ class LSTEventSource(EventSource): EventSource for LST r0 data. """ - + def __init__(self, config=None, tool=None, **kwargs): - + """ Constructor Parameters @@ -63,7 +63,7 @@ def __init__(self, config=None, tool=None, **kwargs): self.multi_file = MultiFiles(self.file_list) - + self.camera_config = self.multi_file.camera_config self.log.info("Read {} input files".format(self.multi_file.num_inputs())) @@ -123,6 +123,10 @@ def _generator(self): @staticmethod def is_compatible(file_path): + from .sst1meventsource import is_fits_in_header + if not is_fits_in_header(file_path): + return False + from astropy.io import fits try: # The file contains two tables: @@ -215,8 +219,8 @@ def fill_r0_camera_container_from_zfile(self, container, event): reshaped_waveform = np.array( event.waveform - ).reshape(n_gains, - self.camera_config.num_pixels, + ).reshape(n_gains, + self.camera_config.num_pixels, container.num_samples) # initialize the waveform container to zero @@ -247,7 +251,7 @@ class MultiFiles: """ This class open all the files in file_list and read the events following - the event_id order + the event_id order """ def __init__(self, file_list): @@ -258,9 +262,9 @@ def __init__(self, file_list): self._camera_config = {} self.camera_config = None - + paths = [] - for file_name in file_list: + for file_name in file_list: paths.append(file_name) Provenance().add_input_file(file_name, role='r0.sub.evt') diff --git a/ctapipe/io/nectarcameventsource.py b/ctapipe/io/nectarcameventsource.py index fea98edf4a4..b57078f782b 100644 --- a/ctapipe/io/nectarcameventsource.py +++ b/ctapipe/io/nectarcameventsource.py @@ -23,8 +23,8 @@ class NectarCAMEventSource(EventSource): """ def __init__(self, config=None, tool=None, **kwargs): - - + + """ Constructor Parameters @@ -87,7 +87,7 @@ def _generator(self): tel_descr = TelescopeDescription(optics, camera) tel_descr.optics.tel_subtype = '' # to correct bug in reading - + self.n_camera_pixels = tel_descr.camera.n_pixels tels = {tel_id: tel_descr} @@ -118,6 +118,10 @@ def _generator(self): @staticmethod def is_compatible(file_path): + from .sst1meventsource import is_fits_in_header + if not is_fits_in_header(file_path): + return False + from astropy.io import fits try: # The file contains two tables: @@ -211,8 +215,8 @@ def fill_r0_camera_container_from_zfile(self, container, event): reshaped_waveform = np.array( event.waveform - ).reshape(n_gains, - self.camera_config.num_pixels, + ).reshape(n_gains, + self.camera_config.num_pixels, container.num_samples) # initialize the waveform container to zero diff --git a/ctapipe/io/sst1meventsource.py b/ctapipe/io/sst1meventsource.py index 35c8a7a99de..e2d80cc5a4b 100644 --- a/ctapipe/io/sst1meventsource.py +++ b/ctapipe/io/sst1meventsource.py @@ -4,6 +4,7 @@ Needs protozfits v1.0.2 from github.com/cta-sst-1m/protozfitsreader """ +import gzip import numpy as np from .eventsource import EventSource from .containers import SST1MDataContainer @@ -12,6 +13,24 @@ __all__ = ['SST1MEventSource'] +def is_fits_in_header(file_path): + '''quick check if file is a FITS file + + by looking into the first 1024 bytes and searching for the string "FITS" + typically used in is_compatible + ''' + # read the first 1kB + with open(file_path, 'rb') as f: + marker_bytes = f.read(1024) + + # if file is gzip, read the first 4 bytes with gzip again + if marker_bytes[0] == 0x1f and marker_bytes[1] == 0x8b: + with gzip.open(file_path, 'rb') as f: + marker_bytes = f.read(1024) + + return b'FITS' in marker_bytes + + class SST1MEventSource(EventSource): def __init__(self, config=None, tool=None, **kwargs): @@ -78,6 +97,9 @@ def _generator(self): @staticmethod def is_compatible(file_path): + if not is_fits_in_header(file_path): + return False + from astropy.io import fits try: h = fits.open(file_path)[1].header From cc1567d808f5330200acd819067339c6d1409f3a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Wed, 23 Jan 2019 14:20:22 +0100 Subject: [PATCH 7/7] Major coordinate refactoring, fixes #899, #900, #505, #416 (#896) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Start refactoring coordinates * Update frame attributes for NominalFrame * Remove weird module check * Fix coordinate usage in muon code * Set argv for tool tests, fixes error when using pytest options * Add roudtrip test * Names to pointings and docstrings. * Renaming in tests and coordinates. no abbreviations. * Fix kwargs * More renames * Use QuantityAttributes were appropriate * Use default=0 for focal_length * More coordinate fixing * More coordinate fixes * Fix coordinate trafo example * Make HorizonFrame an alias to altaz * Add test showing it is possible to go from ICRS to camera frame * Add test that shows separation of planar frame is broken * Add coordinate attribute back * Implement nominal and telescopeframe like skyoffsetframe * Fix imports * Add missing rotation attribute, call it x and y again * Allow transforming of AltAz into AltAz without obstime or location * Use UnitSphericalRepresentation, add tests * Use approx * Use correct representation in camera to telescope * Remove outdated examples * Fix __all__ * Update docstrings * Add example to plot stereo showers * Ugly fix for doc warnings * Fix invalid syntax * Fix style issues * Add test that shows transformation into camera is broken * Wrap TelescopeFrame and NominalFrame angles at 180° * Add comment about mapping function * Rename Telescope and Nominal frame attributes to delta_az, delta_alt, fix camera coordinate trafo * Do not set astropy config at import * Remove comment about iers from docs * Update nominal_frame.py * Update telescope_frame.py * Update docstring of camera coordinate * Start reworking coord docs --- ctapipe/coordinates/__init__.py | 23 +- ctapipe/coordinates/angular_frames.py | 464 ------------------ ctapipe/coordinates/camera_frame.py | 134 +++++ ctapipe/coordinates/ground_frames.py | 63 +-- ctapipe/coordinates/horizon_frame.py | 39 ++ ctapipe/coordinates/nominal_frame.py | 105 ++++ ctapipe/coordinates/representation.py | 8 +- ctapipe/coordinates/telescope_frame.py | 102 ++++ ctapipe/coordinates/tests/test_coordinates.py | 105 +++- .../coordinates/tests/test_nominal_frame.py | 77 +++ ctapipe/coordinates/tests/test_roundtrip.py | 28 ++ .../coordinates/tests/test_telescope_frame.py | 79 +++ ctapipe/image/muon/muon_diagnostic_plots.py | 38 +- ctapipe/image/muon/muon_reco_functions.py | 41 +- .../image/muon/tests/test_muon_ring_finder.py | 2 - ctapipe/reco/HillasReconstructor.py | 42 +- ctapipe/reco/ImPACT.py | 32 +- ctapipe/reco/hillas_intersection.py | 36 +- .../reco/tests/test_HillasReconstructor.py | 39 +- ctapipe/tools/tests/test_tools.py | 5 + docs/coordinates/index.rst | 92 +--- examples/coordinate_transformations.py | 63 ++- examples/plot_array_hillas.py | 4 +- examples/plot_showers_in_nominal.py | 87 ++++ 24 files changed, 974 insertions(+), 734 deletions(-) delete mode 100644 ctapipe/coordinates/angular_frames.py create mode 100644 ctapipe/coordinates/camera_frame.py create mode 100644 ctapipe/coordinates/horizon_frame.py create mode 100644 ctapipe/coordinates/nominal_frame.py create mode 100644 ctapipe/coordinates/telescope_frame.py create mode 100644 ctapipe/coordinates/tests/test_nominal_frame.py create mode 100644 ctapipe/coordinates/tests/test_roundtrip.py create mode 100644 ctapipe/coordinates/tests/test_telescope_frame.py create mode 100644 examples/plot_showers_in_nominal.py diff --git a/ctapipe/coordinates/__init__.py b/ctapipe/coordinates/__init__.py index b0c259c9587..6d336beeb20 100644 --- a/ctapipe/coordinates/__init__.py +++ b/ctapipe/coordinates/__init__.py @@ -1,6 +1,19 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst -""" +''' Coordinates. -""" -from .angular_frames import * -from .ground_frames import * +''' +from .horizon_frame import HorizonFrame +from .telescope_frame import TelescopeFrame +from .nominal_frame import NominalFrame +from .ground_frames import GroundFrame, TiltedGroundFrame, project_to_ground +from .camera_frame import CameraFrame + + +# HorizonFrame is not in __all__ because the docs complain otherwise +__all__ = [ + 'TelescopeFrame', + 'CameraFrame', + 'NominalFrame', + 'GroundFrame', + 'TiltedGroundFrame', + 'project_to_ground', +] diff --git a/ctapipe/coordinates/angular_frames.py b/ctapipe/coordinates/angular_frames.py deleted file mode 100644 index a772a816be2..00000000000 --- a/ctapipe/coordinates/angular_frames.py +++ /dev/null @@ -1,464 +0,0 @@ -"""This module defines the important coordinate systems to be used in -reconstruction with the CTA pipeline and the transformations between -this different systems. Frames and transformations are defined using -the astropy.coordinates framework. This module contains transformations -between angular systems. - -For examples on usage see examples/coordinate_transformations.py - -This code is based on the coordinate transformations performed in the -read_hess code - -TODO: - -- Tests Tests Tests! -- Check cartesian system is still accurate for the nominal and - telescope systems (may need a spherical system) -- Benchmark transformation time -- should use `astropy.coordinates.Angle` for all angles here -""" - -import astropy.units as u -import numpy as np -from astropy.coordinates import (BaseCoordinateFrame, - CartesianRepresentation, - UnitSphericalRepresentation, - FunctionTransform, RepresentationMapping) - -try: - # FrameAttribute was renamed Attribute in astropy 2.0 - # TODO: should really use subclasses like QuantityAttribute - from astropy.coordinates import FrameAttribute as Attribute -except ImportError: - from astropy.coordinates import Attribute - -from astropy.coordinates import frame_transform_graph -from numpy import cos, sin, arctan, arctan2, arcsin, sqrt, arccos, tan -from ..coordinates.representation import PlanarRepresentation - -__all__ = [ - 'CameraFrame', - 'TelescopeFrame', - 'NominalFrame', - 'HorizonFrame' -] - - -class CameraFrame(BaseCoordinateFrame): - """Camera coordinate frame. The camera frame is a simple physical - cartesian frame, describing the 2 dimensional position of objects - in the focal plane of the telescope Most Typically this will be - used to describe the positions of the pixels in the focal plane - - Frame attributes: - - * ``focal_length`` - Focal length of the telescope as a unit quantity (usually meters) - * ``rotation`` - Rotation angle of the camera (0 deg in most cases) - """ - default_representation = PlanarRepresentation - focal_length = Attribute(default=None) - rotation = Attribute(default=0 * u.deg) - pointing_direction = Attribute(default=None) - array_direction = Attribute(default=None) - - -class TelescopeFrame(BaseCoordinateFrame): - """Telescope coordinate frame. Cartesian system to describe the - angular offset of a given position in reference to pointing - direction of a given telescope When pointing corrections become - available they should be applied to the transformation between - this frame and the camera frame - - Frame attributes: - - * ``focal_length`` - Focal length of the telescope as a unit quantity (usually meters) - * ``rotation`` - Rotation angle of the camera (0 deg in most cases) - * ``pointing_direction`` - Alt,Az direction of the telescope pointing - - """ - default_representation = PlanarRepresentation - pointing_direction = Attribute(default=None) - - -class NominalFrame(BaseCoordinateFrame): - """Nominal coordinate frame. Cartesian system to describe the angular - offset of a given position in reference to pointing direction of a - nominal array pointing position. In most cases this frame is the - same as the telescope frame, however in the case of divergent - pointing they will differ. Event reconstruction should be - performed in this system - - Frame attributes: - - * ``array_direction`` - Alt,Az direction of the array pointing - * ``pointing_direction`` - Alt,Az direction of the telescope pointing - - """ - default_representation = PlanarRepresentation - pointing_direction = Attribute(default=None) - array_direction = Attribute(default=None) - - -class HorizonFrame(BaseCoordinateFrame): - """Horizon coordinate frame. Spherical system used to describe the direction - of a given position, in terms of the altitude and azimuth of the system. In - practice this is functionally identical as the astropy AltAz system, but this - implementation allows us to pass array pointing information, allowing us to directly - transform to the Horizon Frame from the Camera system. - The Following attributes are carried over from the telescope frame - to allow a direct transformation from the camera frame - - Frame attributes: - - * ``array_direction`` - Alt,Az direction of the array pointing - * ``pointing_direction`` - Alt,Az direction of the telescope pointing - - """ - default_representation = UnitSphericalRepresentation - - frame_specific_representation_info = { - 'spherical': [RepresentationMapping('lon', 'az'), - RepresentationMapping('lat', 'alt')], - } - - frame_specific_representation_info['unitspherical'] = frame_specific_representation_info['spherical'] - - pointing_direction = Attribute(default=None) - array_direction = Attribute(default=None) - - -# Transformations defined below this point - -def altaz_to_offset(obj_azimuth, obj_altitude, azimuth, altitude): - """ - Function to convert a given altitude and azimuth to a cartesian angular - angular offset with regard to a give reference system - (This function is directly lifted from read_hess) - - Parameters - ---------- - obj_azimuth: float - Event azimuth (radians) - obj_altitude: float - Event altitude (radians) - azimuth: float - Reference system azimuth (radians) - altitude: float - Reference system altitude (radians) - - Returns - ------- - x_off,y_off: (float,float) - Offset of the event in the reference system (in radians) - """ - - diff_az = obj_azimuth - azimuth - cosine_obj_alt = cos(obj_altitude) - - xp0 = -cos(diff_az) * cosine_obj_alt - yp0 = sin(diff_az) * cosine_obj_alt - zp0 = sin(obj_altitude) - - sin_sys_alt = sin(altitude) - cos_sys_alt = cos(altitude) - - xp1 = sin_sys_alt * xp0 + cos_sys_alt * zp0 - yp1 = yp0 - zp1 = -cos_sys_alt * xp0 + sin_sys_alt * zp0 - - disp = tan(arccos(zp1)) - alpha = arctan2(yp1, xp1) - - x_off = disp * cos(alpha) - y_off = disp * sin(alpha) - - return x_off, y_off - - -def offset_to_altaz(x_off, y_off, azimuth, altitude): - """Function to convert an angular offset with regard to a give - reference system to an an absolute altitude and azimuth (This - function is directly lifted from read_hess) - - Parameters - ---------- - x_off: float - X offset of the event in the reference system - y_off: float - Y offset of the event in the reference system - azimuth: float - Reference system azimuth (radians) - altitude: float - Reference system altitude (radians) - - Returns - ------- - obj_altitude,obj_azimuth: (float,float) - Absolute altitude and azimuth of the event - """ - - unit = azimuth.unit - - x_off = x_off.to(u.rad).value - y_off = y_off.to(u.rad).value - azimuth = azimuth.to(u.rad).value - altitude = altitude.to(u.rad).value - - offset = sqrt(x_off * x_off + y_off * y_off) - pos = np.where(offset == 0) # find offset 0 positions - - if len(pos[0]) > 0: - if np.isscalar(offset): - offset = 1e-14 - else: - offset[pos] = 1e-14 # add a very small offset to prevent math errors - - atan_off = arctan(offset) - - sin_atan_off = sin(atan_off) - xp1 = x_off * (sin_atan_off / offset) - yp1 = y_off * (sin_atan_off / offset) - zp1 = cos(atan_off) - - sin_obj_alt = sin(altitude) - cos_obj_alt = cos(altitude) - - xp0 = sin_obj_alt * xp1 - cos_obj_alt * zp1 - yp0 = yp1 - zp0 = cos_obj_alt * xp1 + sin_obj_alt * zp1 - - obj_altitude = arcsin(zp0) - obj_azimuth = arctan2(yp0, -xp0) + azimuth - - if len(pos[0]) > 0: - if np.isscalar(obj_altitude): - obj_altitude = altitude - obj_azimuth = azimuth - else: - obj_altitude[pos] = altitude - obj_azimuth[pos] = azimuth - - obj_altitude = obj_altitude * u.rad - obj_azimuth = obj_azimuth * u.rad - - # if obj_azimuth.value < 0.: - # obj_azimuth += 2.*pi - # elif obj_azimuth.value >= (2.*pi ): - # obj_azimuth -= 2.*pi - - return obj_altitude.to(unit), obj_azimuth.to(unit) - - -# Transformation between nominal and AltAz system - - -@frame_transform_graph.transform(FunctionTransform, NominalFrame, HorizonFrame) -def nominal_to_altaz(norm_coord, altaz_coord): - """ - Transformation from nominal system to astropy AltAz system - - Parameters - ---------- - norm_coord: `astropy.coordinates.SkyCoord` - nominal system - altaz_coord: `astropy.coordinates.SkyCoord` - AltAz system - - Returns - ------- - AltAz Coordinates - """ - alt_norm, az_norm = norm_coord.array_direction.alt, norm_coord.array_direction.az - - if type(norm_coord.x.value).__module__ != np.__name__: - x_off = np.zeros(1) - x_off[0] = norm_coord.x.value - x_off = x_off * norm_coord.x.unit - y_off = np.zeros(1) - y_off[0] = norm_coord.y.value - y_off = y_off * norm_coord.y.unit - else: - x_off = norm_coord.x - y_off = norm_coord.y - - altitude, azimuth = offset_to_altaz(x_off, y_off, az_norm, alt_norm) - - representation = UnitSphericalRepresentation(lon=azimuth.to(u.deg), lat=altitude.to(u.deg)) - return altaz_coord.realize_frame(representation) - - -@frame_transform_graph.transform(FunctionTransform, HorizonFrame, NominalFrame) -def altaz_to_nominal(altaz_coord, norm_coord): - """ - Transformation from astropy AltAz system to nominal system - - Parameters - ---------- - altaz_coord: `astropy.coordinates.SkyCoord` - AltAz system - norm_coord: `astropy.coordinates.SkyCoord` - nominal system - - Returns - ------- - nominal Coordinates - """ - alt_norm, az_norm = norm_coord.array_direction.alt, norm_coord.array_direction.az - azimuth = altaz_coord.az - altitude = altaz_coord.alt - x_off, y_off = altaz_to_offset(azimuth, altitude, az_norm, alt_norm) - x_off = x_off * u.rad - y_off = y_off * u.rad - representation = PlanarRepresentation( - x_off.to(u.deg), y_off.to(u.deg)) - - return norm_coord.realize_frame(representation) - - -# Transformation between telescope and nominal frames - - -@frame_transform_graph.transform(FunctionTransform, TelescopeFrame, NominalFrame) -def telescope_to_nominal(tel_coord, norm_frame): - """ - Coordinate transformation from telescope frame to nominal frame - - Parameters - ---------- - tel_coord: `astropy.coordinates.SkyCoord` - TelescopeFrame system - norm_frame: `astropy.coordinates.SkyCoord` - NominalFrame system - - Returns - ------- - NominalFrame coordinates - """ - alt_tel, az_tel = tel_coord.pointing_direction.alt, tel_coord.pointing_direction.az - alt_norm, az_norm = norm_frame.array_direction.alt, norm_frame.array_direction.az - alt_trans, az_trans = offset_to_altaz( - tel_coord.x, tel_coord.y, az_tel, alt_tel) - - x_off, y_off = altaz_to_offset(az_trans, alt_trans, az_norm, alt_norm) - x_off = x_off * u.rad - y_off = y_off * u.rad - - representation = PlanarRepresentation( - x_off.to(tel_coord.x.unit), y_off.to(tel_coord.x.unit)) - - return norm_frame.realize_frame(representation) - - -@frame_transform_graph.transform(FunctionTransform, NominalFrame, TelescopeFrame) -def nominal_to_telescope(norm_coord, tel_frame): - """ - Coordinate transformation from nominal to telescope system - - Parameters - ---------- - norm_coord: `astropy.coordinates.SkyCoord` - NominalFrame system - tel_frame: `astropy.coordinates.SkyCoord` - TelescopeFrame system - - Returns - ------- - TelescopeFrame coordinates - - """ - alt_tel, az_tel = tel_frame.pointing_direction.alt, tel_frame.pointing_direction.az - alt_norm, az_norm = norm_coord.array_direction.alt, norm_coord.array_direction.az - - alt_trans, az_trans = offset_to_altaz( - norm_coord.x, norm_coord.y, az_norm, alt_norm) - x_off, y_off = altaz_to_offset(az_trans, alt_trans, az_tel, alt_tel) - x_off = x_off * u.rad - y_off = y_off * u.rad - - representation = PlanarRepresentation(x_off.to(norm_coord.x.unit), y_off.to(norm_coord.x.unit)) - - return tel_frame.realize_frame(representation) - - -# Transformations between camera frame and telescope frame -@frame_transform_graph.transform(FunctionTransform, CameraFrame, TelescopeFrame) -def camera_to_telescope(camera_coord, telescope_frame): - """ - Transformation between CameraFrame and TelescopeFrame - - Parameters - ---------- - camera_coord: `astropy.coordinates.SkyCoord` - CameraFrame system - telescope_frame: `astropy.coordinates.SkyCoord` - TelescopeFrame system - Returns - ------- - TelescopeFrame coordinate - """ - x_pos = camera_coord.cartesian.x - y_pos = camera_coord.cartesian.y - - rot = camera_coord.rotation - if rot == 0: - x_rotated = x_pos - y_rotated = y_pos - else: - x_rotated = x_pos * cos(rot) - y_pos * sin(rot) - y_rotated = x_pos * sin(rot) + y_pos * cos(rot) - - focal_length = camera_coord.focal_length - - x_rotated = (x_rotated / focal_length) * u.rad - y_rotated = (y_rotated / focal_length) * u.rad - representation = PlanarRepresentation(x_rotated, y_rotated) - - return telescope_frame.realize_frame(representation) - - -@frame_transform_graph.transform(FunctionTransform, TelescopeFrame, CameraFrame) -def telescope_to_camera(telescope_coord, camera_frame): - """ - Transformation between TelescopeFrame and CameraFrame - - Parameters - ---------- - telescope_coord: `astropy.coordinates.SkyCoord` - TelescopeFrame system - camera_frame: `astropy.coordinates.SkyCoord` - CameraFrame system - - Returns - ------- - CameraFrame Coordinates - """ - x_pos = telescope_coord.cartesian.x - y_pos = telescope_coord.cartesian.y - # reverse the rotation applied to get to this system - rot = camera_frame.rotation * -1 - - if rot == 0: # if no rotation applied save a few cycles - x_rotated = x_pos - y_rotated = y_pos - else: # or else rotate all positions around the camera centre - x_rotated = x_pos * cos(rot) - y_pos * sin(rot) - y_rotated = x_pos * sin(rot) + y_pos * cos(rot) - - focal_length = camera_frame.focal_length - # Remove distance units here as we are using small angle approx - x_rotated = x_rotated.to(u.rad) * (focal_length / u.m) - y_rotated = y_rotated.to(u.rad) * (focal_length / u.m) - - representation = CartesianRepresentation( - x_rotated.value * u.m, y_rotated.value * u.m, 0 * u.m) - - return camera_frame.realize_frame(representation) diff --git a/ctapipe/coordinates/camera_frame.py b/ctapipe/coordinates/camera_frame.py new file mode 100644 index 00000000000..ff1eed03aee --- /dev/null +++ b/ctapipe/coordinates/camera_frame.py @@ -0,0 +1,134 @@ +import numpy as np +import astropy.units as u +from astropy.coordinates import ( + BaseCoordinateFrame, + CoordinateAttribute, + QuantityAttribute, + TimeAttribute, + EarthLocationAttribute, + FunctionTransform, + frame_transform_graph, + CartesianRepresentation, + UnitSphericalRepresentation, +) + +from .horizon_frame import HorizonFrame +from .telescope_frame import TelescopeFrame +from .representation import PlanarRepresentation + + +class CameraFrame(BaseCoordinateFrame): + ''' + Camera coordinate frame. + + The camera frame is a 2d cartesian frame, + describing position of objects in the focal plane of the telescope. + + The frame is defined as in H.E.S.S., starting at the horizon, + the telescope is pointed to magnetic north in azimuth and then up to zenith. + + Now, x points north and y points west, so in this orientation, the + camera coordinates line up with the CORSIKA ground coordinate system. + + MAGIC and FACT use a different camera coordinate system: + Standing at the dish, looking at the camera, x points right, y points up. + To transform MAGIC/FACT to ctapipe, do x' = -y, y' = -x. + + Attributes + ---------- + + focal_length : u.Quantity[length] + Focal length of the telescope as a unit quantity (usually meters) + rotation : u.Quantity[angle] + Rotation angle of the camera (0 deg in most cases) + telescope_pointing : SkyCoord[HorizonFrame] + Pointing direction of the telescope as SkyCoord in HorizonFrame + obstime : Time + Observation time + location : EarthLocation + location of the telescope + ''' + default_representation = PlanarRepresentation + + focal_length = QuantityAttribute(default=0, unit=u.m) + rotation = QuantityAttribute(default=0 * u.deg, unit=u.rad) + telescope_pointing = CoordinateAttribute(frame=HorizonFrame, default=None) + + obstime = TimeAttribute(default=None) + location = EarthLocationAttribute(default=None) + + +@frame_transform_graph.transform(FunctionTransform, CameraFrame, TelescopeFrame) +def camera_to_telescope(camera_coord, telescope_frame): + ''' + Transformation between CameraFrame and TelescopeFrame. + Is called when a SkyCoord is transformed from CameraFrame into TelescopeFrame + ''' + x_pos = camera_coord.cartesian.x + y_pos = camera_coord.cartesian.y + + rot = camera_coord.rotation + if rot == 0: # if no rotation applied save a few cycles + x_rotated = x_pos + y_rotated = y_pos + else: + cosrot = np.cos(rot) + sinrot = np.sin(rot) + x_rotated = x_pos * cosrot - y_pos * sinrot + y_rotated = x_pos * sinrot + y_pos * cosrot + + focal_length = camera_coord.focal_length + + # this assumes an equidistant mapping function of the telescope optics + # or a small angle approximation of most other mapping functions + # this could be replaced by actually defining the mapping function + # as an Attribute of `CameraFrame` that maps f(r, focal_length) -> theta, + # where theta is the angle to the optical axis and r is the distance + # to the camera center in the focal plane + delta_alt = u.Quantity((x_rotated / focal_length).to_value(), u.rad) + delta_az = u.Quantity((y_rotated / focal_length).to_value(), u.rad) + + representation = UnitSphericalRepresentation(lat=delta_alt, lon=delta_az) + + return telescope_frame.realize_frame(representation) + + +@frame_transform_graph.transform(FunctionTransform, TelescopeFrame, CameraFrame) +def telescope_to_camera(telescope_coord, camera_frame): + ''' + Transformation between TelescopeFrame and CameraFrame + + Is called when a SkyCoord is transformed from TelescopeFrame into CameraFrame + ''' + x_pos = telescope_coord.delta_alt + y_pos = telescope_coord.delta_az + # reverse the rotation applied to get to this system + rot = -camera_frame.rotation + + if rot.value == 0.0: # if no rotation applied save a few cycles + x_rotated = x_pos + y_rotated = y_pos + else: # or else rotate all positions around the camera centre + cosrot = np.cos(rot) + sinrot = np.sin(rot) + x_rotated = x_pos * cosrot - y_pos * sinrot + y_rotated = x_pos * sinrot + y_pos * cosrot + + focal_length = camera_frame.focal_length + + # this assumes an equidistant mapping function of the telescope optics + # or a small angle approximation of most other mapping functions + # this could be replaced by actually defining the mapping function + # as an Attribute of `CameraFrame` that maps f(theta, focal_length) -> r, + # where theta is the angle to the optical axis and r is the distance + # to the camera center in the focal plane + x_rotated = x_rotated.to_value(u.rad) * focal_length + y_rotated = y_rotated.to_value(u.rad) * focal_length + + representation = CartesianRepresentation( + x_rotated, + y_rotated, + 0 * u.m + ) + + return camera_frame.realize_frame(representation) diff --git a/ctapipe/coordinates/ground_frames.py b/ctapipe/coordinates/ground_frames.py index d775a3a5d27..c80eed02663 100644 --- a/ctapipe/coordinates/ground_frames.py +++ b/ctapipe/coordinates/ground_frames.py @@ -13,29 +13,24 @@ - Tests Tests Tests! """ - import astropy.units as u import numpy as np -from astropy.coordinates import (BaseCoordinateFrame, - CartesianRepresentation, - FunctionTransform) - -try: - # FrameAttribute was renamed Attribute in astropy 2.0 - # TODO: should really use subclasses like QuantityAttribute - from astropy.coordinates import FrameAttribute as Attribute -except ImportError: - from astropy.coordinates import Attribute - - +from astropy.coordinates import ( + BaseCoordinateFrame, + CartesianRepresentation, + FunctionTransform, + CoordinateAttribute, +) from astropy.coordinates import frame_transform_graph from numpy import cos, sin -from ..coordinates.representation import PlanarRepresentation + +from .representation import PlanarRepresentation +from .horizon_frame import HorizonFrame __all__ = [ 'GroundFrame', 'TiltedGroundFrame', - 'project_to_ground' + 'project_to_ground', ] @@ -50,10 +45,6 @@ class GroundFrame(BaseCoordinateFrame): """ default_representation = CartesianRepresentation - # Pointing direction of the tilted system (alt,az), - # could be the telescope pointing direction or the reconstructed shower - # direction - pointing_direction = Attribute(default=None) class TiltedGroundFrame(BaseCoordinateFrame): @@ -64,7 +55,7 @@ class TiltedGroundFrame(BaseCoordinateFrame): reconstruction of the shower core position Frame attributes: - + * ``pointing_direction`` Alt,Az direction of the tilted reference plane @@ -73,9 +64,7 @@ class TiltedGroundFrame(BaseCoordinateFrame): # Pointing direction of the tilted system (alt,az), # could be the telescope pointing direction or the reconstructed shower # direction - pointing_direction = Attribute(default=None) - -# Transformations defined below this point + pointing_direction = CoordinateAttribute(default=None, frame=HorizonFrame) def get_shower_trans_matrix(azimuth, altitude): @@ -119,28 +108,27 @@ def get_shower_trans_matrix(azimuth, altitude): @frame_transform_graph.transform(FunctionTransform, GroundFrame, TiltedGroundFrame) -def ground_to_tilted(ground_coord, tilted_coord): +def ground_to_tilted(ground_coord, tilted_frame): """ Transformation from ground system to tilted ground system Parameters ---------- ground_coord: `astropy.coordinates.SkyCoord` - GroundFrame system - tilted_coord: `astropy.coordinates.SkyCoord` - TiltedGroundFrame system + Coordinate in GroundFrame + tilted_frame: `ctapipe.coordinates.TiltedFrame` + Frame to transform to Returns ------- - TiltedGroundFrame coordinates + SkyCoordinate transformed to `tilted_frame` coordinates """ x_grd = ground_coord.cartesian.x y_grd = ground_coord.cartesian.y z_grd = ground_coord.cartesian.z - altitude, azimuth = tilted_coord.pointing_direction.alt, tilted_coord.pointing_direction.az - altitude = altitude.to(u.rad) - azimuth = azimuth.to(u.rad) + altitude = tilted_frame.pointing_direction.alt.to(u.rad) + azimuth = tilted_frame.pointing_direction.az.to(u.rad) trans = get_shower_trans_matrix(azimuth, altitude) x_tilt = trans[0][0] * x_grd + trans[0][1] * y_grd + trans[0][2] * z_grd @@ -148,12 +136,12 @@ def ground_to_tilted(ground_coord, tilted_coord): representation = PlanarRepresentation(x_tilt, y_tilt) - return tilted_coord.realize_frame(representation) + return tilted_frame.realize_frame(representation) @frame_transform_graph.transform(FunctionTransform, TiltedGroundFrame, GroundFrame) -def tilted_to_ground(tilted_coord, ground_coord): +def tilted_to_ground(tilted_coord, ground_frame): """ Transformation from tilted ground system to ground system @@ -161,7 +149,7 @@ def tilted_to_ground(tilted_coord, ground_coord): ---------- tilted_coord: `astropy.coordinates.SkyCoord` TiltedGroundFrame system - ground_coord: `astropy.coordinates.SkyCoord` + ground_frame: `astropy.coordinates.SkyCoord` GroundFrame system Returns @@ -171,9 +159,8 @@ def tilted_to_ground(tilted_coord, ground_coord): x_tilt = tilted_coord.x y_tilt = tilted_coord.y - altitude, azimuth = tilted_coord.pointing_direction.alt, tilted_coord.pointing_direction.az - altitude = altitude.to(u.rad) - azimuth = azimuth.to(u.rad) + altitude = tilted_coord.pointing_direction.alt.to(u.rad) + azimuth = tilted_coord.pointing_direction.az.to(u.rad) trans = get_shower_trans_matrix(azimuth, altitude) @@ -183,7 +170,7 @@ def tilted_to_ground(tilted_coord, ground_coord): representation = CartesianRepresentation(x_grd, y_grd, z_grd) - grd = ground_coord.realize_frame(representation) + grd = ground_frame.realize_frame(representation) return grd diff --git a/ctapipe/coordinates/horizon_frame.py b/ctapipe/coordinates/horizon_frame.py new file mode 100644 index 00000000000..8ca69967df6 --- /dev/null +++ b/ctapipe/coordinates/horizon_frame.py @@ -0,0 +1,39 @@ +from astropy.coordinates import ( + AltAz, + FunctionTransformWithFiniteDifference, + CIRS, + frame_transform_graph, +) +import warnings + + +HorizonFrame = AltAz + + +# astropy requires all AltAz to have locations +# and obstimes so they can be converted into true skycoords +# Also it is required to transform one AltAz into another one +# This forbids it to use AltAz without setting location and obstime +# here, the astropy behaviour is defined so that it is assumed, +# that if no information about location or obstime is known, both are the same +@frame_transform_graph.transform( + FunctionTransformWithFiniteDifference, + HorizonFrame, + HorizonFrame +) +def altaz_to_altaz(from_coo, to_frame): + # check if coordinates have obstimes defined + obstime = from_coo.obstime + if from_coo.obstime is None: + warnings.warn('Horizontal coordinate has no obstime, assuming same frame') + obstime = to_frame.obstime + + location = from_coo.location + if from_coo.obstime is None: + warnings.warn('Horizontal coordinate has no location, assuming same frame') + location = to_frame.location + + if obstime is None or location is None: + return to_frame.realize_frame(from_coo.spherical) + + return from_coo.transform_to(CIRS(obstime=obstime)).transform_to(to_frame) diff --git a/ctapipe/coordinates/nominal_frame.py b/ctapipe/coordinates/nominal_frame.py new file mode 100644 index 00000000000..8a4fcf2cd23 --- /dev/null +++ b/ctapipe/coordinates/nominal_frame.py @@ -0,0 +1,105 @@ +''' +The code in this module is basically a copy of +http://docs.astropy.org/en/stable/_modules/astropy/coordinates/builtin_frames/skyoffset.html + +We are just not creating a metaclass and a factory but directly building the +corresponding class. +''' +import astropy.units as u +from astropy.coordinates.matrix_utilities import ( + rotation_matrix, + matrix_product, + matrix_transpose, +) +from astropy.coordinates import ( + frame_transform_graph, + FunctionTransform, + DynamicMatrixTransform, + UnitSphericalRepresentation, + BaseCoordinateFrame, + CoordinateAttribute, + TimeAttribute, + EarthLocationAttribute, + RepresentationMapping, + Angle, +) + +from .horizon_frame import HorizonFrame + + +class NominalFrame(BaseCoordinateFrame): + ''' + Nominal coordinate frame. + + A Frame using a UnitSphericalRepresentation. + This is basically the same as a HorizonCoordinate, but the + origin is at an arbitray position in the sky. + This is what astropy calls a SkyOffsetCoordinate + + If the telescopes are in divergent pointing, this Frame can be + used to transform to a common system. + + Attributes + ---------- + + origin: SkyCoord[HorizonFrame] + Origin of this frame as a HorizonCoordinate + obstime: Tiem + Observation time + location: EarthLocation + Location of the telescope + ''' + frame_specific_representation_info = { + UnitSphericalRepresentation: [ + RepresentationMapping('lon', 'delta_az'), + RepresentationMapping('lat', 'delta_alt'), + ] + } + default_representation = UnitSphericalRepresentation + + origin = CoordinateAttribute(default=None, frame=HorizonFrame) + + obstime = TimeAttribute(default=None) + location = EarthLocationAttribute(default=None) + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + # make sure telescope coordinate is in range [-180°, 180°] + if isinstance(self._data, UnitSphericalRepresentation): + self._data.lon.wrap_angle = Angle(180, unit=u.deg) + + +@frame_transform_graph.transform(FunctionTransform, NominalFrame, NominalFrame) +def skyoffset_to_skyoffset(from_telescope_coord, to_telescope_frame): + """Transform between two skyoffset frames.""" + + intermediate_from = from_telescope_coord.transform_to( + from_telescope_coord.origin + ) + intermediate_to = intermediate_from.transform_to( + to_telescope_frame.origin + ) + return intermediate_to.transform_to(to_telescope_frame) + + +@frame_transform_graph.transform(DynamicMatrixTransform, HorizonFrame, NominalFrame) +def reference_to_skyoffset(reference_frame, telescope_frame): + """Convert a reference coordinate to an sky offset frame.""" + + # Define rotation matrices along the position angle vector, and + # relative to the origin. + origin = telescope_frame.origin.spherical + mat1 = rotation_matrix(-origin.lat, 'y') + mat2 = rotation_matrix(origin.lon, 'z') + return matrix_product(mat1, mat2) + + +@frame_transform_graph.transform(DynamicMatrixTransform, NominalFrame, HorizonFrame) +def skyoffset_to_reference(skyoffset_coord, reference_frame): + """Convert an sky offset frame coordinate to the reference frame""" + + # use the forward transform, but just invert it + mat = reference_to_skyoffset(reference_frame, skyoffset_coord) + # transpose is the inverse because mat is a rotation matrix + return matrix_transpose(mat) diff --git a/ctapipe/coordinates/representation.py b/ctapipe/coordinates/representation.py index a9771dd2743..90ebb4a5785 100644 --- a/ctapipe/coordinates/representation.py +++ b/ctapipe/coordinates/representation.py @@ -1,8 +1,10 @@ """ This module defines any reference systems which may be needed in addition """ - -from astropy.coordinates import BaseRepresentation, CartesianRepresentation +from astropy.coordinates import ( + BaseRepresentation, + CartesianRepresentation, +) import astropy.units as u from collections import OrderedDict from numpy import broadcast_arrays @@ -30,7 +32,6 @@ class PlanarRepresentation(BaseRepresentation): def __init__(self, x, y, copy=True, **kwargs): - if x is None or y is None: raise ValueError( 'x and y are required to instantiate CartesianRepresentation' @@ -57,7 +58,6 @@ def __init__(self, x, y, copy=True, **kwargs): self._y = y self._differentials = {} - @property def x(self): """ diff --git a/ctapipe/coordinates/telescope_frame.py b/ctapipe/coordinates/telescope_frame.py new file mode 100644 index 00000000000..f70b996a078 --- /dev/null +++ b/ctapipe/coordinates/telescope_frame.py @@ -0,0 +1,102 @@ +''' +The code in this module is basically a copy of +http://docs.astropy.org/en/stable/_modules/astropy/coordinates/builtin_frames/skyoffset.html + +We are just not creating a metaclass and a factory but directly building the +corresponding class. +''' +import astropy.units as u +from astropy.coordinates.matrix_utilities import ( + rotation_matrix, + matrix_product, + matrix_transpose, +) +from astropy.coordinates import ( + frame_transform_graph, + FunctionTransform, + DynamicMatrixTransform, + BaseCoordinateFrame, + CoordinateAttribute, + TimeAttribute, + EarthLocationAttribute, + RepresentationMapping, + UnitSphericalRepresentation, + Angle +) + +from .horizon_frame import HorizonFrame + + +class TelescopeFrame(BaseCoordinateFrame): + ''' + Telescope coordinate frame. + + A Frame using a UnitSphericalRepresentation. + This is basically the same as a HorizonCoordinate, but the + origin is at the telescope's pointing direction. + This is what astropy calls a SkyOffsetCoordinate + + Attributes + ---------- + + telescope_pointing: SkyCoord[HorizonFrame] + Coordinate of the telescope pointing in HorizonFrame + obstime: Tiem + Observation time + location: EarthLocation + Location of the telescope + ''' + frame_specific_representation_info = { + UnitSphericalRepresentation: [ + RepresentationMapping('lon', 'delta_az'), + RepresentationMapping('lat', 'delta_alt'), + ] + } + default_representation = UnitSphericalRepresentation + + telescope_pointing = CoordinateAttribute(default=None, frame=HorizonFrame) + + obstime = TimeAttribute(default=None) + location = EarthLocationAttribute(default=None) + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + # make sure telescope coordinate is in range [-180°, 180°] + if isinstance(self._data, UnitSphericalRepresentation): + self._data.lon.wrap_angle = Angle(180, unit=u.deg) + + +@frame_transform_graph.transform(FunctionTransform, TelescopeFrame, TelescopeFrame) +def skyoffset_to_skyoffset(from_telescope_coord, to_telescope_frame): + '''Transform between two skyoffset frames.''' + + intermediate_from = from_telescope_coord.transform_to( + from_telescope_coord.telescope_pointing + ) + intermediate_to = intermediate_from.transform_to( + to_telescope_frame.telescope_pointing + ) + return intermediate_to.transform_to(to_telescope_frame) + + +@frame_transform_graph.transform(DynamicMatrixTransform, HorizonFrame, TelescopeFrame) +def reference_to_skyoffset(reference_frame, telescope_frame): + '''Convert a reference coordinate to an sky offset frame.''' + + # Define rotation matrices along the position angle vector, and + # relative to the telescope_pointing. + telescope_pointing = telescope_frame.telescope_pointing.spherical + mat1 = rotation_matrix(-telescope_pointing.lat, 'y') + mat2 = rotation_matrix(telescope_pointing.lon, 'z') + return matrix_product(mat1, mat2) + + +@frame_transform_graph.transform(DynamicMatrixTransform, TelescopeFrame, HorizonFrame) +def skyoffset_to_reference(skyoffset_coord, reference_frame): + '''Convert an sky offset frame coordinate to the reference frame''' + + # use the forward transform, but just invert it + mat = reference_to_skyoffset(reference_frame, skyoffset_coord) + # transpose is the inverse because mat is a rotation matrix + return matrix_transpose(mat) diff --git a/ctapipe/coordinates/tests/test_coordinates.py b/ctapipe/coordinates/tests/test_coordinates.py index 5ebe4acdaac..cf81e5d79ba 100644 --- a/ctapipe/coordinates/tests/test_coordinates.py +++ b/ctapipe/coordinates/tests/test_coordinates.py @@ -1,6 +1,10 @@ import numpy as np import astropy.units as u -from astropy.coordinates import SkyCoord +from astropy.coordinates import SkyCoord, EarthLocation +from astropy.time import Time +from pytest import approx + +location = EarthLocation.of_site('Roque de los Muchachos') def test_cam_to_nominal(): @@ -9,13 +13,89 @@ def test_cam_to_nominal(): telescope_pointing = SkyCoord(alt=70 * u.deg, az=0 * u.deg, frame=HorizonFrame()) array_pointing = SkyCoord(alt=72 * u.deg, az=0 * u.deg, frame=HorizonFrame()) - cam_frame = CameraFrame(focal_length=28 * u.m, pointing_direction=telescope_pointing) + cam_frame = CameraFrame(focal_length=28 * u.m, telescope_pointing=telescope_pointing) cam = SkyCoord(x=0.5 * u.m, y=0.1 * u.m, frame=cam_frame) - nom_frame = NominalFrame(array_direction=array_pointing) + nom_frame = NominalFrame(origin=array_pointing) cam.transform_to(nom_frame) +def test_icrs_to_camera(): + from ctapipe.coordinates import CameraFrame, HorizonFrame + + obstime = Time('2013-11-01T03:00') + location = EarthLocation.of_site('Roque de los Muchachos') + horizon_frame = HorizonFrame(location=location, obstime=obstime) + + # simulate crab "on" observations + crab = SkyCoord(ra='05h34m31.94s', dec='22d00m52.2s') + telescope_pointing = crab.transform_to(horizon_frame) + + camera_frame = CameraFrame( + focal_length=28 * u.m, + telescope_pointing=telescope_pointing, + location=location, obstime=obstime, + ) + + ceta_tauri = SkyCoord(ra='5h37m38.6854231s', dec='21d08m33.158804s') + ceta_tauri_camera = ceta_tauri.transform_to(camera_frame) + + camera_center = SkyCoord(0 * u.m, 0 * u.m, frame=camera_frame) + crab_camera = crab.transform_to(camera_frame) + + assert crab_camera.x.to_value(u.m) == approx(0.0, abs=1e-10) + assert crab_camera.y.to_value(u.m) == approx(0.0, abs=1e-10) + + # assert ceta tauri is in FoV + assert camera_center.separation_3d(ceta_tauri_camera) < u.Quantity(0.6, u.m) + + +def test_telescope_separation(): + from ctapipe.coordinates import TelescopeFrame, HorizonFrame + + telescope_pointing = SkyCoord( + alt=70 * u.deg, + az=0 * u.deg, + frame=HorizonFrame() + ) + + telescope_frame = TelescopeFrame(telescope_pointing=telescope_pointing) + tel1 = SkyCoord( + delta_az=0 * u.deg, delta_alt=0 * u.deg, frame=telescope_frame + ) + tel2 = SkyCoord( + delta_az=0 * u.deg, delta_alt=1 * u.deg, frame=telescope_frame + ) + + assert tel1.separation(tel2) == u.Quantity(1, u.deg) + + +def test_separation_is_the_same(): + from ctapipe.coordinates import TelescopeFrame, HorizonFrame + + obstime = Time('2013-11-01T03:00') + location = EarthLocation.of_site('Roque de los Muchachos') + horizon_frame = HorizonFrame(location=location, obstime=obstime) + + crab = SkyCoord(ra='05h34m31.94s', dec='22d00m52.2s') + ceta_tauri = SkyCoord(ra='5h37m38.6854231s', dec='21d08m33.158804s') + + # simulate crab "on" observations + telescope_pointing = crab.transform_to(horizon_frame) + + telescope_frame = TelescopeFrame( + telescope_pointing=telescope_pointing, + location=location, + obstime=obstime, + ) + + ceta_tauri_telescope = ceta_tauri.transform_to(telescope_frame) + crab_telescope = crab.transform_to(telescope_frame) + + sep = ceta_tauri_telescope.separation(crab_telescope).to_value(u.deg) + assert ceta_tauri.separation(crab).to_value(u.deg) == approx(sep, rel=1e-4) + + def test_cam_to_tel(): from ctapipe.coordinates import CameraFrame, TelescopeFrame @@ -26,24 +106,23 @@ def test_cam_to_tel(): focal_length = 15 * u.m + camera_frame = CameraFrame(focal_length=focal_length) # first define the camera frame - camera_coord = CameraFrame(pix_x, pix_y, focal_length=focal_length) + camera_coord = SkyCoord(pix_x, pix_y, frame=camera_frame) # then use transform to function to convert to a new system # making sure to give the required values for the conversion # (these are not checked yet) telescope_coord = camera_coord.transform_to(TelescopeFrame()) - assert telescope_coord.x[0] == (1 / 15) * u.rad + assert telescope_coord.delta_az[0] == (1 / 15) * u.rad # check rotation - camera_coord = CameraFrame(pix_x, pix_y, focal_length=focal_length) + camera_coord = SkyCoord(pix_x, pix_y, frame=camera_frame) telescope_coord_rot = camera_coord.transform_to(TelescopeFrame()) - assert telescope_coord_rot.y[0] - (1 / 15) * u.rad < 1e-6 * u.rad + assert telescope_coord_rot.delta_alt[0] - (1 / 15) * u.rad < 1e-6 * u.rad # The Transform back - camera_coord2 = telescope_coord.transform_to( - CameraFrame(focal_length=focal_length) - ) + camera_coord2 = telescope_coord.transform_to(camera_frame) # Check separation assert camera_coord.separation_3d(camera_coord2)[0] == 0 * u.m @@ -54,7 +133,7 @@ def test_ground_to_tilt(): # define ground coordinate grd_coord = GroundFrame(x=1 * u.m, y=2 * u.m, z=0 * u.m) - pointing_direction = HorizonFrame(alt=90 * u.deg, az=0 * u.deg) + pointing_direction = SkyCoord(alt=90 * u.deg, az=0 * u.deg, frame=HorizonFrame()) # Convert to tilted frame at zenith (should be the same) tilt_coord = grd_coord.transform_to( @@ -63,14 +142,14 @@ def test_ground_to_tilt(): assert tilt_coord.separation_3d(grd_coord) == 0 * u.m # Check 180 degree rotation reverses y coordinate - pointing_direction = HorizonFrame(alt=90 * u.deg, az=180 * u.deg) + pointing_direction = SkyCoord(alt=90 * u.deg, az=180 * u.deg, frame=HorizonFrame()) tilt_coord = grd_coord.transform_to( TiltedGroundFrame(pointing_direction=pointing_direction) ) assert np.abs(tilt_coord.y + 2. * u.m) < 1e-5 * u.m # Check that if we look at horizon the x coordinate is 0 - pointing_direction = HorizonFrame(alt=0 * u.deg, az=0 * u.deg) + pointing_direction = SkyCoord(alt=0 * u.deg, az=0 * u.deg, frame=HorizonFrame()) tilt_coord = grd_coord.transform_to( TiltedGroundFrame(pointing_direction=pointing_direction) ) diff --git a/ctapipe/coordinates/tests/test_nominal_frame.py b/ctapipe/coordinates/tests/test_nominal_frame.py new file mode 100644 index 00000000000..d8538593451 --- /dev/null +++ b/ctapipe/coordinates/tests/test_nominal_frame.py @@ -0,0 +1,77 @@ +from astropy.coordinates import SkyCoord +import astropy.units as u +from pytest import approx + + +def test_nominal_to_horizontal_alt0_az0(): + from ctapipe.coordinates.nominal_frame import NominalFrame + from ctapipe.coordinates.horizon_frame import HorizonFrame + + horizon_frame = HorizonFrame() + pointing = SkyCoord(az=0 * u.deg, alt=0 * u.deg, frame=horizon_frame) + + nominal_frame = NominalFrame(origin=pointing) + + nominal_coord = SkyCoord( + delta_az=1 * u.deg, delta_alt=0 * u.deg, frame=nominal_frame + ) + horizon_coord = nominal_coord.transform_to(horizon_frame) + assert horizon_coord.az.deg == 1.0 + assert horizon_coord.alt.deg == 0.0 + + nominal_coord = SkyCoord( + delta_az=-1 * u.deg, delta_alt=0 * u.deg, frame=nominal_frame + ) + horizon_coord = nominal_coord.transform_to(horizon_frame) + assert horizon_coord.az.wrap_at('180d').deg == -1.0 + assert horizon_coord.alt.deg == 0.0 + + nominal_coord = SkyCoord(delta_az=0 * u.deg, delta_alt=1 * u.deg, frame=nominal_frame) + horizon_coord = nominal_coord.transform_to(horizon_frame) + assert horizon_coord.az.deg == 0.0 + assert horizon_coord.alt.deg == 1.0 + + nominal_coord = SkyCoord( + delta_az=0 * u.deg, delta_alt=-1 * u.deg, frame=nominal_frame + ) + horizon_coord = nominal_coord.transform_to(horizon_frame) + assert horizon_coord.az.deg == 0.0 + assert horizon_coord.alt.deg == -1.0 + + +def test_nominal_to_horizontal_alt0_az180(): + from ctapipe.coordinates.nominal_frame import NominalFrame + from ctapipe.coordinates.horizon_frame import HorizonFrame + + horizon_frame = HorizonFrame() + pointing = SkyCoord(az=180 * u.deg, alt=0 * u.deg, frame=horizon_frame) + + nominal_frame = NominalFrame(origin=pointing) + + nominal_coord = SkyCoord( + delta_az=1 * u.deg, delta_alt=0 * u.deg, frame=nominal_frame + ) + horizon_coord = nominal_coord.transform_to(horizon_frame) + assert horizon_coord.az.deg == approx(181.0) + assert horizon_coord.alt.deg == 0.0 + + nominal_coord = SkyCoord( + delta_az=-1 * u.deg, delta_alt=0 * u.deg, frame=nominal_frame + ) + horizon_coord = nominal_coord.transform_to(horizon_frame) + assert horizon_coord.az.deg == approx(179.0) + assert horizon_coord.alt.deg == 0.0 + + nominal_coord = SkyCoord( + delta_az=0 * u.deg, delta_alt=1 * u.deg, frame=nominal_frame + ) + horizon_coord = nominal_coord.transform_to(horizon_frame) + assert horizon_coord.az.deg == 180.0 + assert horizon_coord.alt.deg == 1.0 + + nominal_coord = SkyCoord( + delta_az=0 * u.deg, delta_alt=-1 * u.deg, frame=nominal_frame + ) + horizon_coord = nominal_coord.transform_to(horizon_frame) + assert horizon_coord.az.deg == 180.0 + assert horizon_coord.alt.deg == -1.0 diff --git a/ctapipe/coordinates/tests/test_roundtrip.py b/ctapipe/coordinates/tests/test_roundtrip.py new file mode 100644 index 00000000000..e4b3fa0cefa --- /dev/null +++ b/ctapipe/coordinates/tests/test_roundtrip.py @@ -0,0 +1,28 @@ +from astropy.coordinates import SkyCoord +from pytest import approx +import astropy.units as u + + +def test_roundtrip_camera_horizon(): + from ctapipe.coordinates import CameraFrame, TelescopeFrame, HorizonFrame + + telescope_pointing = SkyCoord(alt=70 * u.deg, az=0 * u.deg, frame=HorizonFrame()) + camera_frame = CameraFrame( + focal_length=28 * u.m, + telescope_pointing=telescope_pointing + ) + + cam_coord = SkyCoord(x=0.5 * u.m, y=0.1 * u.m, frame=camera_frame) + telescope_coord = cam_coord.transform_to(TelescopeFrame()) + horizon_coord = telescope_coord.transform_to(HorizonFrame()) + + back_telescope_coord = horizon_coord.transform_to(TelescopeFrame()) + back_cam_coord = back_telescope_coord.transform_to(camera_frame) + + delta_az = back_telescope_coord.delta_az.to_value(u.deg) + delta_alt = back_telescope_coord.delta_alt.to_value(u.deg) + assert delta_az == approx(telescope_coord.delta_az.to_value(u.deg)) + assert delta_alt == approx(telescope_coord.delta_alt.to_value(u.deg)) + + assert back_cam_coord.x.to_value(u.m) == approx(cam_coord.x.to_value(u.m)) + assert back_cam_coord.y.to_value(u.m) == approx(cam_coord.y.to_value(u.m)) diff --git a/ctapipe/coordinates/tests/test_telescope_frame.py b/ctapipe/coordinates/tests/test_telescope_frame.py new file mode 100644 index 00000000000..885768fe9ed --- /dev/null +++ b/ctapipe/coordinates/tests/test_telescope_frame.py @@ -0,0 +1,79 @@ +from astropy.coordinates import SkyCoord +import astropy.units as u +from pytest import approx + + +def test_telescope_to_horizontal_alt0_az0(): + from ctapipe.coordinates.telescope_frame import TelescopeFrame + from ctapipe.coordinates.horizon_frame import HorizonFrame + + horizon_frame = HorizonFrame() + pointing = SkyCoord(az=0 * u.deg, alt=0 * u.deg, frame=horizon_frame) + + telescope_frame = TelescopeFrame(telescope_pointing=pointing) + + telescope_coord = SkyCoord( + delta_az=1 * u.deg, delta_alt=0 * u.deg, frame=telescope_frame + ) + horizon_coord = telescope_coord.transform_to(horizon_frame) + assert horizon_coord.az.deg == 1.0 + assert horizon_coord.alt.deg == 0.0 + + telescope_coord = SkyCoord( + delta_az=-1 * u.deg, delta_alt=0 * u.deg, frame=telescope_frame + ) + horizon_coord = telescope_coord.transform_to(horizon_frame) + assert horizon_coord.az.wrap_at('180d').deg == -1.0 + assert horizon_coord.alt.deg == 0.0 + + telescope_coord = SkyCoord( + delta_az=0 * u.deg, delta_alt=1 * u.deg, frame=telescope_frame + ) + horizon_coord = telescope_coord.transform_to(horizon_frame) + assert horizon_coord.az.deg == 0.0 + assert horizon_coord.alt.deg == 1.0 + + telescope_coord = SkyCoord( + delta_az=0 * u.deg, delta_alt=-1 * u.deg, frame=telescope_frame + ) + horizon_coord = telescope_coord.transform_to(horizon_frame) + assert horizon_coord.az.deg == 0.0 + assert horizon_coord.alt.deg == -1.0 + + +def test_telescope_to_horizontal_alt0_az180(): + from ctapipe.coordinates.telescope_frame import TelescopeFrame + from ctapipe.coordinates.horizon_frame import HorizonFrame + + horizon_frame = HorizonFrame() + pointing = SkyCoord(az=180 * u.deg, alt=0 * u.deg, frame=horizon_frame) + + telescope_frame = TelescopeFrame(telescope_pointing=pointing) + + telescope_coord = SkyCoord( + delta_az=1 * u.deg, delta_alt=0 * u.deg, frame=telescope_frame + ) + horizon_coord = telescope_coord.transform_to(horizon_frame) + assert horizon_coord.az.deg == approx(181.0) + assert horizon_coord.alt.deg == 0.0 + + telescope_coord = SkyCoord( + delta_az=-1 * u.deg, delta_alt=0 * u.deg, frame=telescope_frame + ) + horizon_coord = telescope_coord.transform_to(horizon_frame) + assert horizon_coord.az.deg == approx(179.0) + assert horizon_coord.alt.deg == 0.0 + + telescope_coord = SkyCoord( + delta_az=0 * u.deg, delta_alt=1 * u.deg, frame=telescope_frame + ) + horizon_coord = telescope_coord.transform_to(horizon_frame) + assert horizon_coord.az.deg == 180.0 + assert horizon_coord.alt.deg == 1.0 + + telescope_coord = SkyCoord( + delta_az=0 * u.deg, delta_alt=-1 * u.deg, frame=telescope_frame + ) + horizon_coord = telescope_coord.transform_to(horizon_frame) + assert horizon_coord.az.deg == 180.0 + assert horizon_coord.alt.deg == -1.0 diff --git a/ctapipe/image/muon/muon_diagnostic_plots.py b/ctapipe/image/muon/muon_diagnostic_plots.py index edab8196b2c..4da234bff2a 100644 --- a/ctapipe/image/muon/muon_diagnostic_plots.py +++ b/ctapipe/image/muon/muon_diagnostic_plots.py @@ -1,17 +1,18 @@ """ -Set of diagnostic plots relating to muons +Set of diagnostic plots relating to muons For generic use with all muon algorithms """ import matplotlib.mlab as mlab import matplotlib.pyplot as plt import numpy as np +from astropy.coordinates import SkyCoord from astropy import units as u from astropy.table import Table from matplotlib import colors from scipy.stats import norm -from ctapipe.coordinates import (CameraFrame, NominalFrame, HorizonFrame) +from ctapipe.coordinates import CameraFrame, NominalFrame, HorizonFrame from ctapipe.image.cleaning import tailcuts_clean from ctapipe.plotting.camera import CameraPlotter from ctapipe.utils.fitshistogram import Histogram @@ -135,11 +136,11 @@ def plot_muon_event(event, muonparams): # Convert to camera frame (centre & radius) altaz = HorizonFrame(alt=event.mc.alt, az=event.mc.az) - ring_nominal = NominalFrame( - x=muonparams['MuonRingParams'][idx].ring_center_x, - y=muonparams['MuonRingParams'][idx].ring_center_y, - array_direction=altaz, - pointing_direction=altaz) + ring_nominal = SkyCoord( + delta_az=muonparams['MuonRingParams'][idx].ring_center_x, + delta_alt=muonparams['MuonRingParams'][idx].ring_center_y, + frame=NominalFrame(origin=altaz) + ) flen = subarray.tel[tel_id].optics.equivalent_focal_length ring_camcoord = ring_nominal.transform_to(CameraFrame( @@ -149,23 +150,26 @@ def plot_muon_event(event, muonparams): centroid = (ring_camcoord.x.value, ring_camcoord.y.value) - ringrad_camcoord = muonparams['MuonRingParams'][idx].ring_radius.to( - u.rad) \ - * flen * 2. # But not FC? + radius = muonparams['MuonRingParams'][idx].ring_radius + ringrad_camcoord = 2 * radius.to(u.rad) * flen # But not FC? px = subarray.tel[tel_id].camera.pix_x py = subarray.tel[tel_id].camera.pix_y - camera_coord = CameraFrame(x=px, y=py, - focal_length=flen, - rotation=geom.pix_rotation) + camera_coord = SkyCoord( + x=px, + y=py, + frame=CameraFrame( + focal_length=flen, + rotation=geom.pix_rotation, + ) + ) nom_coord = camera_coord.transform_to( - NominalFrame(array_direction=altaz, - pointing_direction=altaz) + NominalFrame(origin=altaz) ) - px = nom_coord.x.to(u.deg) - py = nom_coord.y.to(u.deg) + px = nom_coord.delta_az.to(u.deg) + py = nom_coord.delta_alt.to(u.deg) dist = np.sqrt(np.power(px - muonparams['MuonRingParams'][idx].ring_center_x, 2) + np.power(py - muonparams['MuonRingParams'][idx]. ring_center_y, 2)) diff --git a/ctapipe/image/muon/muon_reco_functions.py b/ctapipe/image/muon/muon_reco_functions.py index 6eb0df48ed9..43a08bcfc99 100644 --- a/ctapipe/image/muon/muon_reco_functions.py +++ b/ctapipe/image/muon/muon_reco_functions.py @@ -1,9 +1,10 @@ import logging +import warnings import numpy as np from astropy import log from astropy import units as u -from astropy.coordinates import Angle +from astropy.coordinates import Angle, SkyCoord from astropy.utils.decorators import deprecated from ctapipe.coordinates import CameraFrame, NominalFrame, HorizonFrame @@ -20,7 +21,7 @@ def analyze_muon_event(event): """ - Generic muon event analyzer. + Generic muon event analyzer. Parameters ---------- @@ -29,7 +30,7 @@ def analyze_muon_event(event): Returns ------- - muonringparam, muonintensityparam : MuonRingParameter + muonringparam, muonintensityparam : MuonRingParameter and MuonIntensityParameter container event """ @@ -91,24 +92,32 @@ def analyze_muon_event(event): clean_mask = tailcuts_clean(geom, image, picture_thresh=tailcuts[0], boundary_thresh=tailcuts[1]) - camera_coord = CameraFrame( - x=x, y=y, - focal_length=teldes.optics.equivalent_focal_length, - rotation=geom.pix_rotation - ) # TODO: correct this hack for values over 90 altval = event.mcheader.run_array_direction[1] - if altval > Angle(90*u.deg): - altval = Angle(90*u.deg) + if altval > Angle(90, unit=u.deg): + warnings.warn('Altitude over 90 degrees') + altval = Angle(90, unit=u.deg) + + telescope_pointing = SkyCoord( + alt=altval, + az=event.mcheader.run_array_direction[0], + frame=HorizonFrame() + ) + camera_coord = SkyCoord( + x=x, y=y, + frame=CameraFrame( + focal_length=teldes.optics.equivalent_focal_length, + rotation=geom.pix_rotation, + telescope_pointing=telescope_pointing, + ) + ) - altaz = HorizonFrame(alt=altval, - az=event.mcheader.run_array_direction[0]) nom_coord = camera_coord.transform_to( - NominalFrame(array_direction=altaz, pointing_direction=altaz) + NominalFrame(origin=telescope_pointing) ) - x = nom_coord.x.to(u.deg) - y = nom_coord.y.to(u.deg) + x = nom_coord.delta_az.to(u.deg) + y = nom_coord.delta_alt.to(u.deg) if(cleaning): img = image * clean_mask @@ -270,7 +279,7 @@ def analyze_muon_source(source): log.info("[FUNCTION] {}".format(__name__)) if geom_dict is None: - geom_dict = {} + geom_dict = {} numev = 0 for event in source: # Put a limit on number of events numev += 1 diff --git a/ctapipe/image/muon/tests/test_muon_ring_finder.py b/ctapipe/image/muon/tests/test_muon_ring_finder.py index bca0f1fc05c..09b9de08649 100644 --- a/ctapipe/image/muon/tests/test_muon_ring_finder.py +++ b/ctapipe/image/muon/tests/test_muon_ring_finder.py @@ -72,8 +72,6 @@ def test_ChaudhuriKunduRingFitter(): clean_toy_mask = tailcuts_clean(geom, toymodel_image, boundary_thresh=5, picture_thresh=10) - # camera_coord = CameraFrame(x=x,y=y,z=np.zeros(x.shape)*u.m, - # focal_length = event.inst.optical_foclen[telid], rotation=geom.pix_rotation) muonring = muon_ring_finder.ChaudhuriKunduRingFitter(None) x = np.rad2deg((geom.pix_x.value / 15.) * u.rad) # .value diff --git a/ctapipe/reco/HillasReconstructor.py b/ctapipe/reco/HillasReconstructor.py index 1c301f6cb8c..2503659a114 100644 --- a/ctapipe/reco/HillasReconstructor.py +++ b/ctapipe/reco/HillasReconstructor.py @@ -7,9 +7,11 @@ from ctapipe.reco.reco_algorithms import Reconstructor from ctapipe.io.containers import ReconstructedShowerContainer +from itertools import combinations + from ctapipe.coordinates import HorizonFrame, CameraFrame, GroundFrame, TiltedGroundFrame, project_to_ground from astropy.coordinates import SkyCoord, spherical_to_cartesian, cartesian_to_spherical -from itertools import combinations +import warnings import numpy as np @@ -137,11 +139,11 @@ class are set to np.nan alt = u.Quantity(list(pointing_alt.values())) az = u.Quantity(list(pointing_az.values())) if np.any(alt != alt[0]) or np.any(az != az[0]): - raise ValueError('Divergent pointing not supported') + warnings.warn('Divergent pointing not supported') - pointing_direction = SkyCoord(alt=alt[0], az=az[0], frame='altaz') + telescope_pointing = SkyCoord(alt=alt[0], az=az[0], frame=HorizonFrame()) # core position estimate using a geometric approach - core_pos = self.estimate_core_position(hillas_dict, pointing_direction) + core_pos = self.estimate_core_position(hillas_dict, telescope_pointing) # container class for reconstructed showers result = ReconstructedShowerContainer() @@ -195,6 +197,7 @@ def initialize_hillas_planes( """ self.hillas_planes = {} + horizon_frame = HorizonFrame() for tel_id, moments in hillas_dict.items(): # we just need any point on the main shower axis a bit away from the cog p2_x = moments.x + 0.1 * u.m * np.cos(moments.psi) @@ -204,24 +207,23 @@ def initialize_hillas_planes( pointing = SkyCoord( alt=pointing_alt[tel_id], az=pointing_az[tel_id], - frame='altaz' + frame=horizon_frame, ) - hf = HorizonFrame( - array_direction=pointing, - pointing_direction=pointing - ) - cf = CameraFrame( + camera_frame = CameraFrame( focal_length=focal_length, - array_direction=pointing, - pointing_direction=pointing, + telescope_pointing=pointing ) - cog_coord = SkyCoord(x=moments.x, y=moments.y, frame=cf) - cog_coord = cog_coord.transform_to(hf) + cog_coord = SkyCoord( + x=moments.x, + y=moments.y, + frame=camera_frame, + ) + cog_coord = cog_coord.transform_to(horizon_frame) - p2_coord = SkyCoord(x=p2_x, y=p2_y, frame=cf) - p2_coord = p2_coord.transform_to(hf) + p2_coord = SkyCoord(x=p2_x, y=p2_y, frame=camera_frame) + p2_coord = p2_coord.transform_to(horizon_frame) circle = HillasPlane( p1=cog_coord, @@ -269,7 +271,7 @@ def estimate_direction(self): return result, err_est_dir - def estimate_core_position(self, hillas_dict, pointing_direction): + def estimate_core_position(self, hillas_dict, telescope_pointing): ''' Estimate the core position by intersection the major ellipse lines of each telescope. @@ -277,7 +279,7 @@ def estimate_core_position(self, hillas_dict, pointing_direction): ----------- hillas_dict: dict[HillasContainer] dictionary of hillas moments - pointing_direction: SkyCoord[AltAz] + telescope_pointing: SkyCoord[HorizonFrame] Pointing direction of the array Returns @@ -292,8 +294,8 @@ def estimate_core_position(self, hillas_dict, pointing_direction): z = np.zeros(len(psi)) uvw_vectors = np.column_stack([np.cos(psi).value, np.sin(psi).value, z]) - tilted_frame = TiltedGroundFrame(pointing_direction=pointing_direction) - ground_frame = GroundFrame(pointing_direction=pointing_direction) + tilted_frame = TiltedGroundFrame(pointing_direction=telescope_pointing) + ground_frame = GroundFrame() positions = [ ( diff --git a/ctapipe/reco/ImPACT.py b/ctapipe/reco/ImPACT.py index ba5fb71028b..13dbceb6ddf 100644 --- a/ctapipe/reco/ImPACT.py +++ b/ctapipe/reco/ImPACT.py @@ -7,6 +7,7 @@ import numpy as np import numpy.ma as ma from astropy import units as u +from astropy.coordinates import SkyCoord from iminuit import Minuit from scipy.optimize import minimize, least_squares from scipy.stats import norm @@ -630,12 +631,15 @@ def predict(self, shower_seed, energy_seed): Reconstructed ImPACT shower geometry and energy """ - horizon_seed = HorizonFrame(az=shower_seed.az, alt=shower_seed.alt) - nominal_seed = horizon_seed.transform_to(NominalFrame( - array_direction=self.array_direction)) + horizon_seed = SkyCoord( + az=shower_seed.az, alt=shower_seed.alt, frame=HorizonFrame() + ) + nominal_seed = horizon_seed.transform_to( + NominalFrame(origin=self.array_direction) + ) - source_x = nominal_seed.x.to(u.rad).value - source_y = nominal_seed.y.to(u.rad).value + source_x = nominal_seed.delta_az.to_value(u.rad) + source_y = nominal_seed.delta_alt.to_value(u.rad) ground = GroundFrame(x=shower_seed.core_x, y=shower_seed.core_y, z=0 * u.m) tilted = ground.transform_to( @@ -668,15 +672,19 @@ def predict(self, shower_seed, energy_seed): # Convert the best fits direction and core to Horizon and ground systems and # copy to the shower container - nominal = NominalFrame(x=fit_params[0] * u.rad, - y=fit_params[1] * u.rad, - array_direction=self.array_direction) + nominal = SkyCoord( + x=fit_params[0] * u.rad, + y=fit_params[1] * u.rad, + frame=NominalFrame(origin=self.array_direction) + ) horizon = nominal.transform_to(HorizonFrame()) shower_result.alt, shower_result.az = horizon.alt, horizon.az - tilted = TiltedGroundFrame(x=fit_params[2] * u.m, - y=fit_params[3] * u.m, - pointing_direction=self.array_direction) + tilted = TiltedGroundFrame( + x=fit_params[2] * u.m, + y=fit_params[3] * u.m, + pointing_direction=self.array_direction + ) ground = project_to_ground(tilted) shower_result.core_x = ground.x @@ -913,4 +921,4 @@ def create_seed(source_x, source_y, tilt_x, tilt_y, energy): [0.5, 2] ] - return seed, step, limits \ No newline at end of file + return seed, step, limits diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py index 3ff234037e4..2890e8a17b9 100644 --- a/ctapipe/reco/hillas_intersection.py +++ b/ctapipe/reco/hillas_intersection.py @@ -13,9 +13,11 @@ import astropy.units as u from ctapipe.reco.reco_algorithms import Reconstructor from ctapipe.io.containers import ReconstructedShowerContainer +from ctapipe.instrument import get_atmosphere_profile_functions + +from astropy.coordinates import SkyCoord from ctapipe.coordinates import NominalFrame, HorizonFrame from ctapipe.coordinates import TiltedGroundFrame, project_to_ground -from ctapipe.instrument import get_atmosphere_profile_functions __all__ = [ 'HillasIntersection' @@ -77,32 +79,40 @@ def predict(self, hillas_parameters, tel_x, tel_y, array_direction): err_x *= u.rad err_y *= u.rad - nom = NominalFrame(x=src_x * u.rad, y=src_y * u.rad, - array_direction=array_direction) + nom = SkyCoord( + x=src_x * u.rad, + y=src_y * u.rad, + frame=NominalFrame(array_direction=array_direction) + ) horiz = nom.transform_to(HorizonFrame()) + result = ReconstructedShowerContainer() result.alt, result.az = horiz.alt, horiz.az - tilt = TiltedGroundFrame(x=core_x * u.m, y=core_y * u.m, - pointing_direction=array_direction) + tilt = SkyCoord( + x=core_x * u.m, + y=core_y * u.m, + frame=TiltedGroundFrame(pointing_direction=array_direction), + ) grd = project_to_ground(tilt) result.core_x = grd.x result.core_y = grd.y - x_max = self.reconstruct_xmax(nom.x, nom.y, - tilt.x, tilt.y, - hillas_parameters, - tel_x, tel_y, - 90 * u.deg - array_direction.alt) + x_max = self.reconstruct_xmax( + nom.x, nom.y, + tilt.x, tilt.y, + hillas_parameters, + tel_x, tel_y, + 90 * u.deg - array_direction.alt, + ) - result.core_uncert = np.sqrt(core_err_x * core_err_x - + core_err_y * core_err_y) * u.m + result.core_uncert = np.sqrt(core_err_x**2 + core_err_y**2) * u.m result.tel_ids = [h for h in hillas_parameters.keys()] result.average_size = np.mean([h.intensity for h in hillas_parameters.values()]) result.is_valid = True - src_error = np.sqrt(err_x * err_x + err_y * err_y) + src_error = np.sqrt(err_x**2 + err_y**2) result.alt_uncert = src_error.to(u.deg) result.az_uncert = src_error.to(u.deg) result.h_max = x_max diff --git a/ctapipe/reco/tests/test_HillasReconstructor.py b/ctapipe/reco/tests/test_HillasReconstructor.py index fbfd9b09bf2..6bd4ff7a5f1 100644 --- a/ctapipe/reco/tests/test_HillasReconstructor.py +++ b/ctapipe/reco/tests/test_HillasReconstructor.py @@ -7,6 +7,7 @@ from ctapipe.reco.HillasReconstructor import HillasReconstructor, HillasPlane from ctapipe.utils import get_dataset_path from astropy.coordinates import SkyCoord +from ctapipe.coordinates import HorizonFrame def test_estimator_results(): @@ -14,21 +15,22 @@ def test_estimator_results(): creating some planes pointing in different directions (two north-south, two east-west) and that have a slight position errors (+- 0.1 m in one of the four cardinal directions """ + horizon_frame = HorizonFrame() - p1 = SkyCoord(alt=43 * u.deg, az=45 * u.deg, frame='altaz') - p2 = SkyCoord(alt=47 * u.deg, az=45 * u.deg, frame='altaz') + p1 = SkyCoord(alt=43 * u.deg, az=45 * u.deg, frame=horizon_frame) + p2 = SkyCoord(alt=47 * u.deg, az=45 * u.deg, frame=horizon_frame) circle1 = HillasPlane(p1=p1, p2=p2, telescope_position=[0, 1, 0] * u.m) - p1 = SkyCoord(alt=44 * u.deg, az=90 * u.deg, frame='altaz') - p2 = SkyCoord(alt=46 * u.deg, az=90 * u.deg, frame='altaz') + p1 = SkyCoord(alt=44 * u.deg, az=90 * u.deg, frame=horizon_frame) + p2 = SkyCoord(alt=46 * u.deg, az=90 * u.deg, frame=horizon_frame) circle2 = HillasPlane(p1=p1, p2=p2, telescope_position=[1, 0, 0] * u.m) - p1 = SkyCoord(alt=44.5 * u.deg, az=45 * u.deg, frame='altaz') - p2 = SkyCoord(alt=46.5 * u.deg, az=45 * u.deg, frame='altaz') + p1 = SkyCoord(alt=44.5 * u.deg, az=45 * u.deg, frame=horizon_frame) + p2 = SkyCoord(alt=46.5 * u.deg, az=45 * u.deg, frame=horizon_frame) circle3 = HillasPlane(p1=p1, p2=p2, telescope_position=[0, -1, 0] * u.m) - p1 = SkyCoord(alt=43.5 * u.deg, az=90 * u.deg, frame='altaz') - p2 = SkyCoord(alt=45.5 * u.deg, az=90 * u.deg, frame='altaz') + p1 = SkyCoord(alt=43.5 * u.deg, az=90 * u.deg, frame=horizon_frame) + p2 = SkyCoord(alt=45.5 * u.deg, az=90 * u.deg, frame=horizon_frame) circle4 = HillasPlane(p1=p1, p2=p2, telescope_position=[-1, 0, 0] * u.m) # creating the fit class and setting the the great circle member @@ -47,40 +49,38 @@ def test_h_max_results(): creating some planes pointing in different directions (two north-south, two east-west) and that have a slight position errors (+- 0.1 m in one of the four cardinal directions """ + horizon_frame = HorizonFrame() - p1 = SkyCoord(alt=0 * u.deg, az=45 * u.deg, frame='altaz') - p2 = SkyCoord(alt=0 * u.deg, az=45 * u.deg, frame='altaz') + p1 = SkyCoord(alt=0 * u.deg, az=45 * u.deg, frame=horizon_frame) + p2 = SkyCoord(alt=0 * u.deg, az=45 * u.deg, frame=horizon_frame) circle1 = HillasPlane(p1=p1, p2=p2, telescope_position=[0, 1, 0] * u.m) - p1 = SkyCoord(alt=0 * u.deg, az=90 * u.deg, frame='altaz') - p2 = SkyCoord(alt=0 * u.deg, az=90 * u.deg, frame='altaz') + p1 = SkyCoord(alt=0 * u.deg, az=90 * u.deg, frame=horizon_frame) + p2 = SkyCoord(alt=0 * u.deg, az=90 * u.deg, frame=horizon_frame) circle2 = HillasPlane(p1=p1, p2=p2, telescope_position=[1, 0, 0] * u.m) - p1 = SkyCoord(alt=0 * u.deg, az=45 * u.deg, frame='altaz') - p2 = SkyCoord(alt=0 * u.deg, az=45 * u.deg, frame='altaz') + p1 = SkyCoord(alt=0 * u.deg, az=45 * u.deg, frame=horizon_frame) + p2 = SkyCoord(alt=0 * u.deg, az=45 * u.deg, frame=horizon_frame) circle3 = HillasPlane(p1=p1, p2=p2, telescope_position=[0, -1, 0] * u.m) - p1 = SkyCoord(alt=0 * u.deg, az=90 * u.deg, frame='altaz') - p2 = SkyCoord(alt=0 * u.deg, az=90 * u.deg, frame='altaz') + p1 = SkyCoord(alt=0 * u.deg, az=90 * u.deg, frame=horizon_frame) + p2 = SkyCoord(alt=0 * u.deg, az=90 * u.deg, frame=horizon_frame) circle4 = HillasPlane(p1=p1, p2=p2, telescope_position=[-1, 0, 0] * u.m) # creating the fit class and setting the the great circle member fit = HillasReconstructor() fit.hillas_planes = {1: circle1, 2: circle2, 3: circle3, 4: circle4} - # performing the direction fit with the minimisation algorithm # and a seed that is perpendicular to the up direction h_max_reco = fit.estimate_h_max() print("h max fit test minimise:", h_max_reco) - # the results should be close to the direction straight up np.testing.assert_allclose(h_max_reco.value, 0, atol=1e-8) # np.testing.assert_allclose(fitted_core_position.value, [0, 0], atol=1e-3) - def test_reconstruction(): """ a test of the complete fit procedure on one event including: @@ -133,4 +133,3 @@ def test_reconstruction(): fit_result.az.to(u.deg) fit_result.core_x.to(u.m) assert fit_result.is_valid - return diff --git a/ctapipe/tools/tests/test_tools.py b/ctapipe/tools/tests/test_tools.py index fe0f252799a..98113e70bae 100644 --- a/ctapipe/tools/tests/test_tools.py +++ b/ctapipe/tools/tests/test_tools.py @@ -1,3 +1,4 @@ +import sys from ctapipe.tools.camdemo import CameraDemo from ctapipe.tools.dump_triggers import DumpTriggersTool from ctapipe.tools.dump_instrument import DumpInstrumentTool @@ -15,6 +16,7 @@ def test_info(): def test_dump_triggers(tmpdir): + sys.argv = ['dump_triggers'] outfile = tmpdir.join("triggers.fits") tool = DumpTriggersTool( @@ -28,6 +30,7 @@ def test_dump_triggers(tmpdir): def test_dump_instrument(tmpdir): + sys.argv = ['dump_instrument'] tmpdir.chdir() tool = DumpInstrumentTool( @@ -41,6 +44,7 @@ def test_dump_instrument(tmpdir): def test_camdemo(): + sys.argv = ['camera_demo'] tool = CameraDemo() tool.num_events = 10 tool.cleanframes = 2 @@ -49,6 +53,7 @@ def test_camdemo(): def test_bokeh_file_viewer(): + sys.argv = ['bokeh_file_viewer'] tool = BokehFileViewer(disable_server=True) tool.run() diff --git a/docs/coordinates/index.rst b/docs/coordinates/index.rst index 55886ae819d..9aeb9099529 100644 --- a/docs/coordinates/index.rst +++ b/docs/coordinates/index.rst @@ -11,100 +11,30 @@ Introduction `ctapipe.coordinates` contains coordinate frame definitions and coordinate transformation routines that are associated with the -reconstruction of Cherenkov telescope events. It is built on -`astropy.coordinates`, which internally use the ERFA coordinate -transformation library, the open-source-licensed fork of the IAU SOFA -system. +reconstruction of Cherenkov telescope events. +It is built on `astropy.coordinates`, +which internally use the ERFA coordinate transformation library, +the open-source-licensed fork of the IAU SOFA system. -Maybe we'll use the `gwcs `__ package -for generalised world coordinate transformations to implement some coordinate -transformations, e.g. precision pointing corrections using a pointing model. - -.. note:: - - Due to the underlying implementation of astropy coordinates and time - representations, the transformation of a single coordinate between - frames or time systems has a small overhead in speed. - This can be mitigated by transforming many coordinates at once, which is - much faster. All routines can accept arrays instead of single values. - - -Earth Orientation Corrections -============================== - -Polar motion and leap-second corrections are taken into account automatically -by the use of `astropy.coordinates` and `astropy.time`. Updates to the IERS -tables are loaded automatically if an internet connection is available. - - - -TODO: ------ - -- include time transform examples (using `astropy.time`), particularly - to the CTA-standard Mission-Elapsed (MET) time system, which is a - Terrestrial Time with a specific MJD_offset. Getting Started =============== -The coordinate library defines a set of *frames*, which represent -different coordinate reprentations. Coordinates should be described -using an `astropy.coordinates.SkyCoord` in the appropriate frame. +The coordinate library defines a set of *frames*, +which represent different coordinate reprentations. +Coordinates should be described using an `astropy.coordinates.SkyCoord` in the appropriate frame. The following special frames are defined for CTA: * `CameraFrame` -* `GroundFrame` -* `TiltedGroundFrame` +* `TelescopeFrame` * `NominalFrame` * `HorizonFrame` -* `TelescopeFrame` +* `GroundFrame` +* `TiltedGroundFrame` they can be transformed to and from any other `astropy.coordinates` frame, like -`astropy.coordinates.AltAz` or `astropy.coordinates.ICRS` (RA/Dec) - -Examples --------- - -.. code-block:: python - - >>> import astropy.units as u - >>> from ctapipe.coordinates import CameraFrame, TelescopeFrame - >>> camera_coord = CameraFrame(x=1*u.m, y=2*u.m, z=0*u.m) - >>> print(camera_coord) - - >>> telescope_coord = camera_coord.transform_to(TelescopeFrame) - >>> print(telescope_coord) - - - -Note that Coordinate objects can accept arrays as input, so if one -wants to transform a large list of coordinates from one frame to -another, it is most efficient to pass them all at once. - -.. code-block:: python - - import astropy.units as u - from ctapipe.coordinates import CameraFrame, TelescopeFrame - x = np.array([0.1,0.2,0.3,0.4,0.5]) * u.m - y = np.array([-0.1,0.4,-0.4,0.6,-0.6]) * u.m - camera_coord = CameraFrame(x=x, y=y, z=0*u.m) - telescope_coord = camera_coord.transform_to(TelescopeFrame) - print(telescope_coord) - print("X:",telescope_coord.x) - print("Y:",telescope_coord.y) - -The result is a set of coordinates:: - - - X: [ 0.00666667 0.01333333 0.02 0.02666667 0.03333333] - Y: [-0.00666667 0.02666667 -0.02666667 0.04 -0.04 ] +`astropy.coordinates.ICRS` (RA/Dec) diff --git a/examples/coordinate_transformations.py b/examples/coordinate_transformations.py index 6d1e00f2de1..1350a1f3d7f 100755 --- a/examples/coordinate_transformations.py +++ b/examples/coordinate_transformations.py @@ -5,10 +5,18 @@ """ import astropy.units as u +from astropy.coordinates import SkyCoord import numpy as np -from ctapipe.coordinates import CameraFrame, TelescopeFrame, GroundFrame, \ - TiltedGroundFrame, NominalFrame, HorizonFrame, project_to_ground +from ctapipe.coordinates import ( + CameraFrame, + TelescopeFrame, + GroundFrame, + TiltedGroundFrame, + NominalFrame, + HorizonFrame, + project_to_ground, +) # Convert camera frame to telescope frame @@ -20,8 +28,11 @@ def cam_to_tel(): # e.g. in this case the position on pixels in the camera pix_x = np.ones(2048) * u.m pix_y = np.ones(2048) * u.m + # first define the camera frame - camera_coord = CameraFrame(pix_x, pix_y, focal_length=15 * u.m, rotation=0 * u.deg) + camera_frame = CameraFrame(focal_length=15 * u.m) + # create a coordinate in that frame + camera_coord = SkyCoord(pix_x, pix_y, frame=camera_frame) # then use transform to function to convert to a new system making sure # to give the required values for the conversion (these are not checked @@ -32,9 +43,7 @@ def cam_to_tel(): print("Telescope Coordinate", telescope_coord) # Transforming back is then easy - camera_coord2 = telescope_coord.transform_to( - CameraFrame(focal_length=15 * u.m, rotation=0 * u.deg) - ) + camera_coord2 = telescope_coord.transform_to(camera_frame) # We can easily check the distance between 2 coordinates in the same frame # In this case they should be the same @@ -46,37 +55,37 @@ def cam_to_tel(): def cam_to_nom(): pix_x = np.ones(2048) * u.m pix_y = np.ones(2048) * u.m - camera_coord = CameraFrame(pix_x, pix_y, focal_length=15 * u.m) - # In this case we bypass the telescope system - nom_coord = camera_coord.transform_to( - NominalFrame( - pointing_direction=HorizonFrame(alt=70 * u.deg, az=180 * u.deg), - array_direction=HorizonFrame(alt=75 * u.deg, az=180 * u.deg) - ) + + pointing_direction = SkyCoord(alt=70 * u.deg, az=180 * u.deg, frame=HorizonFrame()) + camera_frame = CameraFrame( + focal_length=15 * u.m, + telescope_pointing=pointing_direction ) - alt_az = camera_coord.transform_to( - HorizonFrame( - pointing_direction=HorizonFrame(alt=70 * u.deg, az=180 * u.deg), - array_direction=HorizonFrame(alt=75 * u.deg, az=180 * u.deg) - ) + camera_coord = SkyCoord(pix_x, pix_y, frame=camera_frame) + + # In this case we bypass the telescope system + nominal_frame = NominalFrame( + origin=HorizonFrame(alt=75 * u.deg, az=180 * u.deg) ) + nom_coord = camera_coord.transform_to(nominal_frame) + + horizon = camera_coord.transform_to(HorizonFrame()) print("Nominal Coordinate", nom_coord) - print("AltAz coordinate", alt_az) + print("Horizon coordinate", horizon) # Once we are at the nominal system where most reconstruction will be done we # can then convert to AltAz (currently we cannot transform directly from camera) def nominal_to_altaz(): - t = np.zeros(10) - t[5] = 1 - nom = NominalFrame( - x=t * u.deg, - y=t * u.deg, - array_direction=HorizonFrame(alt=75 * u.deg, az=180 * u.deg) + + nom = SkyCoord( + x=0 * u.deg, + y=0 * u.deg, + frame=NominalFrame(origin=HorizonFrame(alt=75 * u.deg, az=180 * u.deg)) ) - alt_az = nom.transform_to(HorizonFrame) - print("AltAz Coordinate", alt_az) + alt_az = nom.transform_to(HorizonFrame()) + print("HorizonCoordinate", alt_az) # Provided we know when and where the AltAz was measured we can them # convert this to any astronomical diff --git a/examples/plot_array_hillas.py b/examples/plot_array_hillas.py index 2583edfe2c5..56b39c2a213 100644 --- a/examples/plot_array_hillas.py +++ b/examples/plot_array_hillas.py @@ -10,7 +10,7 @@ from astropy import units as u from ctapipe.calib import CameraCalibrator -from ctapipe.coordinates import TiltedGroundFrame +from ctapipe.coordinates import TiltedGroundFrame, HorizonFrame from ctapipe.image import hillas_parameters, tailcuts_clean, \ HillasParameterizationError from ctapipe.io import event_source @@ -61,7 +61,7 @@ # GroundFrame) point_dir = SkyCoord( *event.mcheader.run_array_direction, - frame='altaz' + frame=HorizonFrame() ) tiltedframe = TiltedGroundFrame(pointing_direction=point_dir) if markers: diff --git a/examples/plot_showers_in_nominal.py b/examples/plot_showers_in_nominal.py new file mode 100644 index 00000000000..49e636837e5 --- /dev/null +++ b/examples/plot_showers_in_nominal.py @@ -0,0 +1,87 @@ +from astropy.coordinates import SkyCoord +import astropy.units as u +from ctapipe.io import event_source +from ctapipe.image.cleaning import tailcuts_clean +from ctapipe.calib import CameraCalibrator +from ctapipe.utils.datasets import get_dataset_path +import matplotlib.pyplot as plt +import numpy as np + +from ctapipe.coordinates import HorizonFrame, CameraFrame, NominalFrame + + +cleaning_level = { + 'LSTCam': (3.5, 7.5, 2), # ?? (3, 6) for Abelardo... + 'FlashCam': (4, 8, 2), # there is some scaling missing? + 'ASTRICam': (5, 7, 2), +} + + +input_url = get_dataset_path('gamma_test_large.simtel.gz') + + +with event_source(input_url=input_url) as source: + calibrator = CameraCalibrator( + eventsource=source, + ) + + for event in source: + + calibrator.calibrate(event) + + nominal_frame = NominalFrame( + origin=SkyCoord(alt=70 * u.deg, az=0 * u.deg, frame=HorizonFrame) + ) + + nom_delta_az = [] + nom_delta_alt = [] + photons = [] + + for tel_id, dl1 in event.dl1.tel.items(): + camera = event.inst.subarray.tels[tel_id].camera + focal_length = event.inst.subarray.tels[tel_id].optics.equivalent_focal_length + image = dl1.image[0] + + # telescope mc info + mc_tel = event.mc.tel[tel_id] + + telescope_pointing = SkyCoord( + alt=mc_tel['altitude_raw'], + az=mc_tel['azimuth_raw'], + unit='rad', frame=HorizonFrame(), + ) + camera_frame = CameraFrame( + telescope_pointing=telescope_pointing, focal_length=focal_length + ) + + boundary, picture, min_neighbors = cleaning_level[camera.cam_id] + clean = tailcuts_clean( + camera, + image, + boundary_thresh=boundary, + picture_thresh=picture, + min_number_picture_neighbors=min_neighbors + ) + + cam_coords = SkyCoord( + camera.pix_x[clean], + camera.pix_y[clean], + frame=camera_frame + ) + nom = cam_coords.transform_to(nominal_frame) + nom_delta_az.append(nom.delta_az.to_value(u.deg)) + nom_delta_alt.append(nom.delta_alt.to_value(u.deg)) + photons.append(image[clean]) + + nom_delta_az = np.concatenate(nom_delta_az) + nom_delta_alt = np.concatenate(nom_delta_alt) + photons = np.concatenate(photons) + + nom_delta_az = np.repeat(nom_delta_az, photons.astype(int)) + nom_delta_alt = np.repeat(nom_delta_alt, photons.astype(int)) + + plt.hexbin(nom_delta_az, nom_delta_alt, gridsize=50, extent=[-5, 5, -5, 5]) + plt.xlabel('delta_az / deg') + plt.ylabel('delta_alt / deg') + plt.gca().set_aspect(1) + plt.show()