From 54d2e3c48aa9912a96cf8041fa889dffe9a505b0 Mon Sep 17 00:00:00 2001 From: vuillaut Date: Thu, 9 Jul 2020 20:27:49 +0200 Subject: [PATCH 1/7] add functions to read instrument info from DL1 file --- lstchain/io/io.py | 237 +++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 236 insertions(+), 1 deletion(-) diff --git a/lstchain/io/io.py b/lstchain/io/io.py index 396db6aea8..9d4fc82c57 100644 --- a/lstchain/io/io.py +++ b/lstchain/io/io.py @@ -16,6 +16,8 @@ from tqdm import tqdm from ctapipe.tools.stage1 import Stage1ProcessorTool from astropy.utils import deprecated +from ctapipe.instrument import OpticsDescription, CameraGeometry, CameraDescription, CameraReadout, \ + TelescopeDescription, SubarrayDescription __all__ = ['read_simu_info_hdf5', @@ -420,7 +422,7 @@ def write_array_info(subarray, output_filename): serialize_meta=serialize_meta, ) - +@deprecated('09/07/2020', message='will be removed in lstchain v0.7') def read_array_info(filename): """ Read array information from HDF5 file. @@ -443,6 +445,239 @@ def read_array_info(filename): return array_info +def read_optic(filename, telescope_name): + """ + Read a specific telescope optics from a DL1 file + + Parameters + ---------- + filename: str + telescope_name: str + + Returns + ------- + `ctapipe.instrument.optics.OpticsDescription` + """ + from astropy.units import Quantity + telescope_optics_path = "/configuration/instrument/telescope/optics" + telescope_optic_table = Table.read(filename, path=telescope_optics_path) + row = telescope_optic_table[np.where(telescope_name == telescope_optic_table['name'])[0][0]] + optics_description = OpticsDescription( + name=row['name'], + num_mirrors=row['num_mirrors'], + equivalent_focal_length=row['equivalent_focal_length'] * telescope_optic_table['equivalent_focal_length'].unit, + mirror_area=row['mirror_area'] * telescope_optic_table['mirror_area'].unit, + num_mirror_tiles=Quantity(row['num_mirror_tiles']), + ) + return optics_description + + +def read_optics(filename): + """ + Read all telescope optics from a DL1 file + + Parameters + ---------- + filename: str + + Returns + ------- + dictionnary of ctapipe.instrument.optics.OpticsDescription by telescope names + """ + telescope_optics_path = "/configuration/instrument/telescope/optics" + telescope_optics_table = Table.read(filename, path=telescope_optics_path) + optics_dict = {} + for telescope_name in telescope_optics_table['name']: + optics_dict[telescope_name] = read_optic(filename, telescope_name) + return optics_dict + + +def read_camera_geometry(filename, camera_name): + """ + Read a specific camera geometry from a DL1 file + + Parameters + ---------- + filename: str + camera_name: str + + Returns + ------- + `ctapipe.instrument.camera.geometry.CameraGeometry` + """ + camera_geometry_path = f"/configuration/instrument/telescope/camera/geometry_{camera_name}" + camera_geometry = CameraGeometry.from_table(Table.read(filename, camera_geometry_path)) + return camera_geometry + + +def read_camera_geometries(filename): + """ + Read all camera geometries from a DL1 file + + Parameters + ---------- + filename + + Returns + ------- + dictionnary of `ctapipe.instrument.camera.geometry.CameraGeometry` by camera name + """ + subarray_layout_path = 'configuration/instrument/subarray/layout' + camera_geoms = {} + for camera_name in set(Table.read(filename, path=subarray_layout_path)['camera_type']): + camera_geoms[camera_name] = read_camera_geometry(filename, camera_name) + return camera_geoms + + +def read_camera_readout(filename, camera_name): + """ + Read a specific camera readout from a DL1 file + + Parameters + ---------- + filename: str + camera_name: str + + Returns + ------- + `ctapipe.instrument.camera.readout.CameraReadout` + """ + camera_readout_path = f"/configuration/instrument/telescope/camera/readout_{camera_name}" + return CameraReadout.from_table(Table.read(filename, path=camera_readout_path)) + + +def read_camera_readouts(filename): + """ + Read all camera readouts from a DL1 file + + Parameters + ---------- + filename: str + + Returns + ------- + dict of `ctapipe.instrument.camera.description.CameraDescription` by tel_id + """ + subarray_layout_path = 'configuration/instrument/subarray/layout' + camera_readouts = {} + for row in Table.read(filename, path=subarray_layout_path): + camera_name = row['camera_type'] + camera_readouts[row['tel_id']] = read_camera_readout(filename, camera_name) + return camera_readouts + + +def read_camera_description(filename, camera_name): + """ + Read a specific camera description from a DL1 file + + Parameters + ---------- + filename: str + camera_name: str + + Returns + ------- + `ctapipe.instrument.camera.description.CameraDescription` + """ + geom = read_camera_geometry(filename, camera_name) + readout = read_camera_readout(filename, camera_name) + return CameraDescription(camera_name, geometry=geom, readout=readout) + + +def read_telescope_description(filename, telescope_name, telescope_type, camera_name): + """ + Read a specific telescope description from a DL1 file + + Parameters + ---------- + filename: str + telescope_name: str + camera_name: str + + Returns + ------- + `ctapipe.instrument.telescope.TelescopeDescription` + """ + optics = read_optic(filename, telescope_name) + camera_descr = read_camera_description(filename, camera_name) + return TelescopeDescription(telescope_name, telescope_type, optics=optics, camera=camera_descr) + + +def read_subarray_table(filename): + """ + Read the subarray as a table from a DL1 file + + Parameters + ---------- + filename: str + + Returns + ------- + `astropy.table.table.Table` + """ + subarray_layout_path = 'configuration/instrument/subarray/layout' + return Table.read(filename, path=subarray_layout_path) + + +def read_telescopes_descriptions(filename): + """ + Read telescopes descriptions from DL1 file + + Parameters + ---------- + filename: str + + Returns + ------- + dict of `ctapipe.instrument.telescope.TelescopeDescription` by tel_id + """ + subarray_table = read_subarray_table(filename) + descriptions = {} + for row in subarray_table: + tel_name = row['name'] + camera_type = row['camera_type'] + optics = read_optic(filename, tel_name) + camera = read_camera_description(filename, camera_type) + descriptions[row['tel_id']] = TelescopeDescription(row['name'], row['type'], optics=optics, camera=camera) + return descriptions + + +def read_telescopes_positions(filename): + """ + Read telescopes positions from DL1 file + + Parameters + ---------- + filename: str + + Returns + ------- + dictionnary of telescopes positions by tel_id + """ + subarray_table = read_subarray_table(filename) + pos_dict = {} + for row in subarray_table['tel_id', 'pos_x', 'pos_y', 'pos_z'].iterrows(): + pos_dict[row[0]] = row[1], row[2], row[3] + return pos_dict + + +def read_subarray_description(filename, subarray_name='LST-1'): + """ + Read subarray description from an HDF5 DL1 file + + Parameters + ---------- + filename: str + + Returns + ------- + `ctapipe.instrument.subarray.SubarrayDescription` + """ + tel_pos = read_telescopes_positions(filename) + tel_descrp = read_telescopes_descriptions(filename) + return SubarrayDescription(subarray_name, tel_positions=tel_pos, tel_descriptions=tel_descrp) + + def check_mcheader(mcheader1, mcheader2): """ Check that the information in two mcheaders are physically consistent. From dcf82f04b5fe802ef24a2ae760ff18464a5b052a Mon Sep 17 00:00:00 2001 From: vuillaut Date: Thu, 9 Jul 2020 21:31:39 +0200 Subject: [PATCH 2/7] tel pos with unit --- lstchain/io/io.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/lstchain/io/io.py b/lstchain/io/io.py index 9d4fc82c57..0e30c12752 100644 --- a/lstchain/io/io.py +++ b/lstchain/io/io.py @@ -656,8 +656,9 @@ def read_telescopes_positions(filename): """ subarray_table = read_subarray_table(filename) pos_dict = {} - for row in subarray_table['tel_id', 'pos_x', 'pos_y', 'pos_z'].iterrows(): - pos_dict[row[0]] = row[1], row[2], row[3] + pos_unit = subarray_table['pos_x'].unit + for row in subarray_table: + pos_dict[row['tel_id']] = np.array([row['pos_x'], row['pos_y'], row['pos_z']]) * pos_unit return pos_dict From 53b41cae245ba6d97d56746f5a47bf6762ebb11a Mon Sep 17 00:00:00 2001 From: vuillaut Date: Thu, 9 Jul 2020 21:32:52 +0200 Subject: [PATCH 3/7] add unit test for subarray info --- lstchain/io/tests/test_io.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/lstchain/io/tests/test_io.py b/lstchain/io/tests/test_io.py index 3361930d57..b485c8ab19 100644 --- a/lstchain/io/tests/test_io.py +++ b/lstchain/io/tests/test_io.py @@ -6,7 +6,7 @@ import os from astropy.table import Table -from lstchain.tests.test_lstchain import dl1_file, test_dir +from lstchain.tests.test_lstchain import mc_gamma_testfile, dl1_file, test_dir merged_dl1_file = os.path.join(test_dir, 'dl1_merged.h5') @@ -106,3 +106,12 @@ def test_trigger_type_in_dl1_params(): from lstchain.io.io import dl1_params_lstcam_key params = pd.read_hdf(dl1_file, key=dl1_params_lstcam_key) assert 'trigger_type' in params.columns + + +@pytest.mark.run(after='test_r0_to_dl1') +def test_read_subarray_description(): + from lstchain.io.io import read_subarray_description + from ctapipe.io import event_source + source = event_source(mc_gamma_testfile) + dl1_subarray = read_subarray_description(dl1_file) + assert(dl1_subarray.info() == source.subarray.info()) \ No newline at end of file From 0979be346d8e319d4aeff6ee917135fecf9686fd Mon Sep 17 00:00:00 2001 From: vuillaut Date: Thu, 9 Jul 2020 21:33:39 +0200 Subject: [PATCH 4/7] add peek in unit test --- lstchain/io/tests/test_io.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lstchain/io/tests/test_io.py b/lstchain/io/tests/test_io.py index b485c8ab19..6219dd0125 100644 --- a/lstchain/io/tests/test_io.py +++ b/lstchain/io/tests/test_io.py @@ -114,4 +114,5 @@ def test_read_subarray_description(): from ctapipe.io import event_source source = event_source(mc_gamma_testfile) dl1_subarray = read_subarray_description(dl1_file) - assert(dl1_subarray.info() == source.subarray.info()) \ No newline at end of file + dl1_subarray.peek() + assert(dl1_subarray.info() == source.subarray.info()) From a6fb58b96429358912c714b15e8050336c8d8025 Mon Sep 17 00:00:00 2001 From: vuillaut Date: Thu, 9 Jul 2020 21:36:14 +0200 Subject: [PATCH 5/7] rename functions with more explicit names (add single) --- lstchain/io/io.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/lstchain/io/io.py b/lstchain/io/io.py index 0e30c12752..0dd52667bb 100644 --- a/lstchain/io/io.py +++ b/lstchain/io/io.py @@ -445,7 +445,7 @@ def read_array_info(filename): return array_info -def read_optic(filename, telescope_name): +def read_single_optics(filename, telescope_name): """ Read a specific telescope optics from a DL1 file @@ -488,11 +488,11 @@ def read_optics(filename): telescope_optics_table = Table.read(filename, path=telescope_optics_path) optics_dict = {} for telescope_name in telescope_optics_table['name']: - optics_dict[telescope_name] = read_optic(filename, telescope_name) + optics_dict[telescope_name] = read_single_optics(filename, telescope_name) return optics_dict -def read_camera_geometry(filename, camera_name): +def read_single_camera_geometry(filename, camera_name): """ Read a specific camera geometry from a DL1 file @@ -525,11 +525,11 @@ def read_camera_geometries(filename): subarray_layout_path = 'configuration/instrument/subarray/layout' camera_geoms = {} for camera_name in set(Table.read(filename, path=subarray_layout_path)['camera_type']): - camera_geoms[camera_name] = read_camera_geometry(filename, camera_name) + camera_geoms[camera_name] = read_single_camera_geometry(filename, camera_name) return camera_geoms -def read_camera_readout(filename, camera_name): +def read_single_camera_readout(filename, camera_name): """ Read a specific camera readout from a DL1 file @@ -562,11 +562,11 @@ def read_camera_readouts(filename): camera_readouts = {} for row in Table.read(filename, path=subarray_layout_path): camera_name = row['camera_type'] - camera_readouts[row['tel_id']] = read_camera_readout(filename, camera_name) + camera_readouts[row['tel_id']] = read_single_camera_readout(filename, camera_name) return camera_readouts -def read_camera_description(filename, camera_name): +def read_single_camera_description(filename, camera_name): """ Read a specific camera description from a DL1 file @@ -579,12 +579,12 @@ def read_camera_description(filename, camera_name): ------- `ctapipe.instrument.camera.description.CameraDescription` """ - geom = read_camera_geometry(filename, camera_name) - readout = read_camera_readout(filename, camera_name) + geom = read_single_camera_geometry(filename, camera_name) + readout = read_single_camera_readout(filename, camera_name) return CameraDescription(camera_name, geometry=geom, readout=readout) -def read_telescope_description(filename, telescope_name, telescope_type, camera_name): +def read_single_telescope_description(filename, telescope_name, telescope_type, camera_name): """ Read a specific telescope description from a DL1 file @@ -598,8 +598,8 @@ def read_telescope_description(filename, telescope_name, telescope_type, camera_ ------- `ctapipe.instrument.telescope.TelescopeDescription` """ - optics = read_optic(filename, telescope_name) - camera_descr = read_camera_description(filename, camera_name) + optics = read_single_optics(filename, telescope_name) + camera_descr = read_single_camera_description(filename, camera_name) return TelescopeDescription(telescope_name, telescope_type, optics=optics, camera=camera_descr) @@ -636,8 +636,8 @@ def read_telescopes_descriptions(filename): for row in subarray_table: tel_name = row['name'] camera_type = row['camera_type'] - optics = read_optic(filename, tel_name) - camera = read_camera_description(filename, camera_type) + optics = read_single_optics(filename, tel_name) + camera = read_single_camera_description(filename, camera_type) descriptions[row['tel_id']] = TelescopeDescription(row['name'], row['type'], optics=optics, camera=camera) return descriptions From f80d8893404cb55871a4b392fa0f4d430d8ebece Mon Sep 17 00:00:00 2001 From: vuillaut Date: Thu, 9 Jul 2020 22:09:36 +0200 Subject: [PATCH 6/7] unit test for subarray --- lstchain/io/tests/test_io.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lstchain/io/tests/test_io.py b/lstchain/io/tests/test_io.py index 6219dd0125..03773d51b6 100644 --- a/lstchain/io/tests/test_io.py +++ b/lstchain/io/tests/test_io.py @@ -115,4 +115,5 @@ def test_read_subarray_description(): source = event_source(mc_gamma_testfile) dl1_subarray = read_subarray_description(dl1_file) dl1_subarray.peek() - assert(dl1_subarray.info() == source.subarray.info()) + dl1_subarray.info() + assert len(dl1_subarray.tels) == len(source.subarray.tels) From 298efef298dea1c87b77dc4136e5956aeb6c8a0a Mon Sep 17 00:00:00 2001 From: vuillaut Date: Thu, 9 Jul 2020 22:21:11 +0200 Subject: [PATCH 7/7] unit test for subarray --- lstchain/io/tests/test_io.py | 1 + 1 file changed, 1 insertion(+) diff --git a/lstchain/io/tests/test_io.py b/lstchain/io/tests/test_io.py index 03773d51b6..f9379654e8 100644 --- a/lstchain/io/tests/test_io.py +++ b/lstchain/io/tests/test_io.py @@ -117,3 +117,4 @@ def test_read_subarray_description(): dl1_subarray.peek() dl1_subarray.info() assert len(dl1_subarray.tels) == len(source.subarray.tels) + assert (dl1_subarray.to_table() == source.subarray.to_table()).all()