diff --git a/.gitignore b/.gitignore index e098b72a3dc..631313313e1 100644 --- a/.gitignore +++ b/.gitignore @@ -12,6 +12,7 @@ __pycache__ ctapipe/_version_cache.py ctapipe/_version.py + # Ignore .c files by default to avoid including generated code. If you want to # add a non-generated .c extension, use `git add -f filename.c`. *.c @@ -25,6 +26,7 @@ MANIFEST # Sphinx docs/api docs/_build +docs/auto_examples # Editors and IDEs @@ -86,7 +88,4 @@ distribute-*.tar.gz target .mypy_cache -examples/notebooks/*.html -examples/notebooks/*.png - provenance.log diff --git a/docs/Makefile b/docs/Makefile index 387d531e5bc..3bf5e372d5f 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -19,7 +19,7 @@ help: @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) clean: - rm -rf api + rm -rf api auto_examples @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/conf.py b/docs/conf.py index 013d575931e..29d5cdee0c5 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -4,7 +4,7 @@ # ctapipe documentation build configuration file, created by # sphinx-quickstart on Fri Jan 6 10:22:58 2017. # -# This file is execfile()d with the current directory set to its +# Thi file is execfile()d with the current directory set to its # containing dir. # # Note that not all possible configuration values are present in this @@ -29,6 +29,9 @@ # Get configuration information from setup.cfg from configparser import ConfigParser +# Sphinx gallery +from sphinx_gallery.sorting import ExplicitOrder, FileNameSortKey + import ctapipe setup_cfg = ConfigParser() @@ -55,6 +58,7 @@ "numpydoc", "sphinx_design", "IPython.sphinxext.ipython_console_highlighting", + "sphinx_gallery.gen_gallery", ] @@ -142,6 +146,31 @@ def setup(app): ("py:class", "ctapipe.compat.StrEnum"), ] +# Sphinx gallery config + +sphinx_gallery_conf = { + "examples_dirs": [ + "../examples", + ], # path to your example scripts + "subsection_order": ExplicitOrder( + [ + "../examples/tutorials", + "../examples/algorithms", + "../examples/core", + "../examples/visualization", + ] + ), + "within_subsection_order": FileNameSortKey, + "nested_sections": False, + "filename_pattern": r".*\.py", + "copyfile_regex": r".*\.png", + "promote_jupyter_magic": True, + "line_numbers": True, + "default_thumb_file": "_static/ctapipe_logo.png", + "pypandoc": True, + "matplotlib_animations": True, +} + # The suffix(es) of source filenames. # You can specify multiple suffix as a list of string: @@ -195,6 +224,12 @@ def setup(app): ".DS_Store", "**.ipynb_checkpoints", "changes", + "user-guide/examples/*/README.rst", + "user-guide/examples/README.rst", + "auto_examples/index.rst", + "auto_examples/*/*.py.md5", + "auto_examples/*/*.py", + "auto_examples/*/*.ipynb", ] # The name of the Pygments (syntax highlighting) style to use. @@ -209,7 +244,7 @@ def setup(app): # Define the json_url for our version switcher. json_url = "https://ctapipe.readthedocs.io/en/latest/_static/switcher.json" -# Define the version we use for matching in the version switcher. +# Define the version we use for matching in the version switcher., version_match = os.getenv("READTHEDOCS_VERSION") # If READTHEDOCS_VERSION doesn't exist, we're not on RTD # If it is an integer, we're in a PR build and the version isn't correct. diff --git a/docs/conftest.py b/docs/conftest.py new file mode 100644 index 00000000000..4c7e9da09ce --- /dev/null +++ b/docs/conftest.py @@ -0,0 +1,5 @@ +collect_ignore = [ + "conf.py", + "_build", + "auto_examples", +] diff --git a/docs/user-guide/examples/InstrumentDescription.ipynb b/docs/user-guide/examples/InstrumentDescription.ipynb deleted file mode 100644 index 2c3b6bd7ea0..00000000000 --- a/docs/user-guide/examples/InstrumentDescription.ipynb +++ /dev/null @@ -1,374 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Working with Instrumental Descriptions\n", - "\n", - "the instrumental description is loaded by the event source, and consists of a hierarchy of classes in the ctapipe.instrument module, the base of which is the `SubarrayDescription`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.utils.datasets import get_dataset_path\n", - "from ctapipe.io import EventSource\n", - "import numpy as np\n", - "\n", - "filename = get_dataset_path(\"gamma_prod5.simtel.zst\")\n", - "\n", - "with EventSource(filename, max_events=1) as source:\n", - " subarray = source.subarray" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## the SubarrayDescription:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.info()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.to_table()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "You can also get a table of just the `OpticsDescriptions` (`CameraGeometry` is more complex and can't be stored on a single table row, so each one can be converted to a table separately)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.to_table(kind=\"optics\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Make a sub-array with only SC-type telescopes:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "sc_tels = [tel_id for tel_id, tel in subarray.tel.items() if tel.optics.n_mirrors == 2]\n", - "newsub = subarray.select_subarray(sc_tels, name=\"SCTels\")\n", - "newsub.info()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "can also do this by using `Table.group_by`" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Explore some of the details of the telescopes" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tel = subarray.tel[1]\n", - "tel" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tel.optics.mirror_area" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tel.optics.n_mirror_tiles" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tel.optics.equivalent_focal_length" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tel.camera" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tel.camera.geometry.pix_x" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib inline\n", - "from ctapipe.visualization import CameraDisplay\n", - "\n", - "CameraDisplay(tel.camera.geometry)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "CameraDisplay(subarray.tel[98].camera.geometry)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Plot the subarray\n", - "\n", - "We'll make a subarray by telescope type and plot each separately, so they appear in different colors. We also calculate the radius using the mirror area (and exagerate it a bit).\n", - "\n", - "This is just for debugging and info, for any \"real\" use, a `visualization.ArrayDisplay` should be used" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.peek()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.footprint" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Get info about the subarray in general" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.telescope_types" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.camera_types" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.optics_types" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from astropy.coordinates import SkyCoord\n", - "from ctapipe.coordinates import GroundFrame\n", - "\n", - "center = SkyCoord(\"10.0 m\", \"2.0 m\", \"0.0 m\", frame=\"groundframe\")\n", - "coords = subarray.tel_coords # a flat list of coordinates by tel_index\n", - "coords.separation(center)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Telescope IDs vs Indices\n", - "\n", - "Note that `subarray.tel` is a dict mapped by `tel_id` (the indentifying number of a telescope). It is possible to have telescope IDs that do not start at 0, are not contiguouous (e.g. if a subarray is selected). Some functions and properties like `tel_coords` are numpy arrays (not dicts) so they are not mapped to the telescope ID, but rather the *index* within this SubarrayDescription. To convert between the two concepts you can do:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.tel_ids_to_indices([1, 5, 23])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "or you can get the indexing array directly in numpy or dict form:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.tel_index_array" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.tel_index_array[[1, 5, 23]]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.tel_indices[\n", - " 1\n", - "] # this is a dict of tel_id -> tel_index, so we can only do one at once" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ids = subarray.get_tel_ids_for_type(subarray.telescope_types[0])\n", - "ids" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "idx = subarray.tel_ids_to_indices(ids)\n", - "idx" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.tel_coords[idx]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "so, with that method you can quickly get many telescope positions at once (the alternative is to use the dict `positions` which maps `tel_id` to a position on the ground" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.positions[1]" - ] - } - ], - "metadata": { - "language_info": { - "codemirror_mode": { - "name": "ipython" - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/docs/user-guide/examples/Tools.ipynb b/docs/user-guide/examples/Tools.ipynb deleted file mode 100644 index 02eb2b8f3c1..00000000000 --- a/docs/user-guide/examples/Tools.ipynb +++ /dev/null @@ -1,618 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Creating command-line Tools" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.core import Tool, Component, TelescopeComponent\n", - "from ctapipe.core.traits import (\n", - " Integer,\n", - " Float,\n", - " List,\n", - " Dict,\n", - " Unicode,\n", - " TraitError,\n", - " observe,\n", - " FloatTelescopeParameter,\n", - " Path,\n", - ")\n", - "import logging\n", - "from time import sleep\n", - "from astropy import units as u\n", - "from ctapipe.utils import get_dataset_path" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "GAMMA_FILE = get_dataset_path(\"gamma_prod5.simtel.zst\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "see https://github.com/ipython/traitlets/blob/master/examples/myapp.py" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setup:\n", - "\n", - "Create a few `Component`s that we will use later in a `Tool`:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "class MyComponent(Component):\n", - " \"\"\"A Component that does stuff\"\"\"\n", - "\n", - " value = Integer(default_value=-1, help=\"Value to use\").tag(config=True)\n", - "\n", - " def do_thing(self):\n", - " self.log.debug(\"Did thing\")\n", - "\n", - "\n", - "# in order to have 2 of the same components at once\n", - "class SecondaryMyComponent(MyComponent):\n", - " \"\"\"A second component\"\"\"\n", - "\n", - " pass\n", - "\n", - "\n", - "class AdvancedComponent(Component):\n", - " \"\"\"An advanced technique\"\"\"\n", - "\n", - " value1 = Integer(default_value=-1, help=\"Value to use\").tag(config=True)\n", - " infile = Path(\n", - " help=\"input file name\",\n", - " exists=None, # set to True to require existing, False for requiring non-existing\n", - " directory_ok=False,\n", - " ).tag(config=True)\n", - " outfile = Path(help=\"output file name\", exists=False, directory_ok=False).tag(\n", - " config=True\n", - " )\n", - "\n", - " def __init__(self, config=None, parent=None, **kwargs):\n", - " super().__init__(config=config, parent=parent, **kwargs)\n", - " # components can have sub components, but these must have\n", - " # then parent=self as argument and be assigned as member\n", - " # so the full config can be received later\n", - " self.subcompent = MyComponent(parent=self)\n", - "\n", - " @observe(\"outfile\")\n", - " def on_outfile_changed(self, change):\n", - " self.log.warning(\"Outfile was changed to '{}'\".format(change))\n", - "\n", - "\n", - "class TelescopeWiseComponent(TelescopeComponent):\n", - " \"\"\"a component that contains parameters that are per-telescope configurable\"\"\"\n", - "\n", - " param = FloatTelescopeParameter(\n", - " help=\"Something configurable with telescope patterns\", default_value=5.0\n", - " ).tag(config=True)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "MyComponent()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "AdvancedComponent(infile=\"test.foo\", outfile=\"out.foo\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "`TelescopeComponents` need to have a subarray given to them in order to work (since they need one to turn a `TelescopeParameter` into a concrete list of values for each telescope. Here we will give a dummy one:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.instrument import SubarrayDescription, TelescopeDescription\n", - "\n", - "subarray = SubarrayDescription.read(GAMMA_FILE)\n", - "subarray.info()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "TelescopeWiseComponent(subarray=subarray)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This TelescopeParameters can then be set using a list of patterns like:\n", - "```python\n", - "component.param = [ \n", - " (\"type\", \"LST*\",3.0), \n", - " (\"type\", \"MST*\", 2.0), \n", - " (id, 25, 4.0) \n", - "]\n", - "```\n", - "\n", - "These get translated into per-telescope-id values once the subarray is registered. After that one acccess the per-telescope id values via:\n", - "```python\n", - "component.param.tel[tel_id]\n", - "```" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Now create an executable Tool that contains the Components\n", - "Note that all the components we wish to be configured via the tool must be added to the `classes` attribute." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "class MyTool(Tool):\n", - " name = \"mytool\"\n", - " description = \"do some things and stuff\"\n", - " aliases = dict(\n", - " infile=\"AdvancedComponent.infile\",\n", - " outfile=\"AdvancedComponent.outfile\",\n", - " iterations=\"MyTool.iterations\",\n", - " )\n", - "\n", - " # Which classes are registered for configuration\n", - " classes = [\n", - " MyComponent,\n", - " AdvancedComponent,\n", - " SecondaryMyComponent,\n", - " TelescopeWiseComponent,\n", - " ]\n", - "\n", - " # local configuration parameters\n", - " iterations = Integer(5, help=\"Number of times to run\", allow_none=False).tag(\n", - " config=True\n", - " )\n", - "\n", - " def setup(self):\n", - " self.comp = MyComponent(parent=self)\n", - " self.comp2 = SecondaryMyComponent(parent=self)\n", - " self.comp3 = TelescopeWiseComponent(parent=self, subarray=subarray)\n", - " self.advanced = AdvancedComponent(parent=self)\n", - "\n", - " def start(self):\n", - " self.log.info(\"Performing {} iterations...\".format(self.iterations))\n", - " for ii in range(self.iterations):\n", - " self.log.info(\"ITERATION {}\".format(ii))\n", - " self.comp.do_thing()\n", - " self.comp2.do_thing()\n", - " sleep(0.1)\n", - "\n", - " def finish(self):\n", - " self.log.warning(\"Shutting down.\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Get Help info\n", - "\n", - "The following allows you to print the help info within a Jupyter notebook, but this same inforamtion would be displayed if the user types:\n", - "```\n", - " mytool --help\n", - "```" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool = MyTool()\n", - "tool" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool.print_help()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The following is equivalant to the user typing `mytool --help-all`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool.print_help(classes=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Run the tool\n", - "\n", - "here we pass in argv since it is a Notebook, but if argv is not specified it's read from `sys.argv`, so the following is the same as running:\n", - "\n", - "```sh\n", - "mytool --log_level=INFO --infile gamma_test.simtel.gz --iterations=3\n", - "```\n", - "\n", - "As Tools are intended to be exectutables, they are raising `SystemExit` on exit.\n", - "Here, we use them to demonstrate how it would work, so we catch the `SystemExit`." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "try:\n", - " tool.run(argv=[\"--infile\", str(GAMMA_FILE), \"--outfile\", \"out.csv\"])\n", - "except SystemExit as e:\n", - " assert e.code == 0, f\"Tool returned with error status {e}\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool.log_format = \"%(asctime)s : %(levelname)s [%(name)s %(funcName)s] %(message)s\"\n", - "\n", - "\n", - "try:\n", - " tool.run(\n", - " argv=[\n", - " \"--log-level\",\n", - " \"INFO\",\n", - " \"--infile\",\n", - " str(GAMMA_FILE),\n", - " \"--outfile\",\n", - " \"out.csv\",\n", - " \"--iterations\",\n", - " \"3\",\n", - " ]\n", - " )\n", - "except SystemExit as e:\n", - " assert e.code == 0, f\"Tool returned with error status {e}\"" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "here we change the log-level to DEBUG:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "try:\n", - " tool.run(\n", - " argv=[\n", - " \"--log-level\",\n", - " \"DEBUG\",\n", - " \"--infile\",\n", - " str(GAMMA_FILE),\n", - " \"--outfile\",\n", - " \"out.csv\",\n", - " ]\n", - " )\n", - "except SystemExit as e:\n", - " assert e.code == 0, f\"Tool returned with error status {e}\"" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "you can also set parameters directly in the class, rather than using the argument/configfile parser. This is useful if you are calling the Tool from a script rather than the command-line" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool.iterations = 1\n", - "tool.log_level = 0\n", - "\n", - "try:\n", - " tool.run([\"--infile\", str(GAMMA_FILE), \"--outfile\", \"out.csv\"])\n", - "except SystemExit as e:\n", - " assert e.code == 0, f\"Tool returned with error status {e}\"" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "see what happens when a value is set that is not of the correct type:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "try:\n", - " tool.iterations = \"badval\"\n", - "except TraitError as E:\n", - " print(\"bad value:\", E)\n", - "except SystemExit as e:\n", - " assert e.code == 0, f\"Tool returned with error status {e}\"" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Example of what happens when you change a parameter that is being \"observed\" in a class. It's handler is called:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool.advanced.outfile = \"Another.txt\"" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "we see that the handler for `outfile` was called, and it receive a change dict that shows the old and new values." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "create a tool using a config file:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool2 = MyTool()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "try:\n", - " tool2.run(argv=[\"--config\", \"Tools.json\"])\n", - "except SystemExit as e:\n", - " assert e.code == 0, f\"Tool returned with error status {e}\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(tool2.advanced.infile)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(tool2.config)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool2.is_setup" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool3 = MyTool()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool3.is_setup" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool3.initialize(argv=[])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool3.is_setup" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool3" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool.setup()\n", - "tool" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool.comp2" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Getting the configuration of an instance" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool.get_current_config()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tool.iterations = 12\n", - "tool.get_current_config()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Writing a Sample Config File" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(tool.generate_config_file())" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.16" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/docs/user-guide/examples/array_display.ipynb b/docs/user-guide/examples/array_display.ipynb deleted file mode 100644 index 6e7ca08f78e..00000000000 --- a/docs/user-guide/examples/array_display.ipynb +++ /dev/null @@ -1,377 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "fb302dc9-65c4-41b0-b8d3-865c0f7845a7", - "metadata": {}, - "source": [ - "# Array Displays\n", - "\n", - "Like `CameraDisplays`, ctapipe provides a way to display information related to the array on the ground: `ArrayDisplay`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "1217a94e-f541-4dfe-8f4f-ec9dd70296fd", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "from ctapipe.visualization import ArrayDisplay\n", - "from ctapipe.instrument import SubarrayDescription\n", - "from ctapipe.coordinates import GroundFrame, EastingNorthingFrame\n", - "from ctapipe.containers import HillasParametersContainer\n", - "\n", - "from astropy.coordinates import SkyCoord\n", - "from astropy import units as u\n", - "\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "\n", - "plt.rcParams[\"figure.figsize\"] = (8, 6)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c7c6775d-ce49-4d83-a3c6-e7f5ed55abfd", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "tel_ids = list(range(1, 5)) + list(range(5, 20)) # just LSTs + one set of MSTs\n", - "\n", - "subarray = SubarrayDescription.read(\n", - " \"dataset://gamma_20deg_0deg_run1___cta-prod5-lapalma_desert-2158m-LaPalma-dark_100evts.simtel.zst\"\n", - ").select_subarray(tel_ids)" - ] - }, - { - "cell_type": "markdown", - "id": "d3b103fc-6b8c-4233-939a-09c5f5d0e60d", - "metadata": {}, - "source": [ - "An array display is created for example in `subarray.peek()`:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a32c4646-9c74-4665-bfaa-81aa1de88dc8", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "subarray.peek()" - ] - }, - { - "cell_type": "markdown", - "id": "a27ed92a-268c-404c-8b96-039a739ea26a", - "metadata": {}, - "source": [ - "However, you can make one manually with a bit more flexibility:" - ] - }, - { - "cell_type": "markdown", - "id": "503dea3f-fff3-4cd2-b938-c0d524a6069c", - "metadata": {}, - "source": [ - "## Constructing an ArrayDisplay" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e27d5c28-dfc7-424b-b511-bcb0a1f35729", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "disp = ArrayDisplay(subarray)" - ] - }, - { - "cell_type": "markdown", - "id": "ea297d1c-4777-4945-81d0-e7807605b742", - "metadata": {}, - "source": [ - "You can specify the Frame you want as long as it is compatible with `GroundFrame`. `EastingNorthingFrame` is probably the most useful. You can also add telescope labels" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "46acf3bd-6cf3-42bb-8a76-163b44d3cc23", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "disp = ArrayDisplay(subarray, frame=EastingNorthingFrame())\n", - "disp.add_labels()" - ] - }, - { - "cell_type": "markdown", - "id": "2925e1f6-eb29-43a7-bda1-c647424e7cf7", - "metadata": {}, - "source": [ - "## Using color to show information \n", - "\n", - "By default the color of the telescope circles correlates to telescope type. However, you can use color to convey other information by setting the `values` attribute, like a trigger pattern" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "141ef1eb-82ca-4e45-a3ae-2b11828c601c", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "plt.set_cmap(\"rainbow\") # the array display will use the current colormap for values\n", - "\n", - "ad = ArrayDisplay(subarray)\n", - "ad.telescopes.set_linewidth(0) # to turn off the telescope borders\n", - "\n", - "trigger_pattern = np.zeros(subarray.n_tels)\n", - "trigger_pattern[\n", - " [\n", - " 1,\n", - " 4,\n", - " 5,\n", - " 6,\n", - " ]\n", - "] = 1\n", - "ad.values = trigger_pattern # display certain telescopes in a color\n", - "ad.add_labels()" - ] - }, - { - "cell_type": "markdown", - "id": "a2c0cec8-46b8-4ea4-b4a5-70d943ed474c", - "metadata": {}, - "source": [ - "or for example, you could use color to represent the telescope distance to the impact point" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8422c9da-df65-44c4-b8bf-da49d95b16fb", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "shower_impact = SkyCoord(200 * u.m, -200 * u.m, 0 * u.m, frame=EastingNorthingFrame())" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e9e87e72-d000-4c14-99eb-665849d83e56", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "plt.set_cmap(\"rainbow\") # the array display will use the current colormap for values\n", - "ad = ArrayDisplay(subarray)\n", - "ad.telescopes.set_linewidth(0) # to turn off the telescope borders\n", - "plt.scatter(shower_impact.easting, shower_impact.northing, marker=\"+\", s=200)\n", - "\n", - "distances = np.hypot(\n", - " subarray.tel_coords.cartesian.x - shower_impact.cartesian.x,\n", - " subarray.tel_coords.cartesian.y - shower_impact.cartesian.y,\n", - ")\n", - "ad.values = distances\n", - "plt.colorbar(ad.telescopes, label=\"Distance (m)\")" - ] - }, - { - "cell_type": "markdown", - "id": "359935df-f551-40cf-8432-b191986c8213", - "metadata": {}, - "source": [ - "## Overlaying vectors\n", - "\n", - "For plotting reconstruction quantities, it's useful to overlay vectors on the telescope positions. `ArrayDisplay` provides functions:\n", - "* `set_vector_uv` to set by cartesian coordinates from the center of each telescope\n", - "* `set_vector_rho_phi` to set by polar coorinates from the center of each telescope\n", - "* `set_vector_hillas` to set vectors from a `dict[int,HillasParameters]` mapping tel_id (not index!) to a set of parameters. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "040a0f4c-ffee-4a25-a419-7d90f5f16a6f", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "np.random.seed(0)\n", - "phis = np.random.uniform(0, 180.0, size=subarray.n_tels) * u.deg\n", - "rhos = np.ones(subarray.n_tels) * 200 * u.m\n", - "\n", - "\n", - "ad = ArrayDisplay(subarray, frame=EastingNorthingFrame(), tel_scale=2)\n", - "ad.set_vector_rho_phi(rho=rhos, phi=phis)" - ] - }, - { - "cell_type": "markdown", - "id": "980be561-28f6-41e6-94cd-ec19d260fe3b", - "metadata": { - "tags": [] - }, - "source": [ - "## Overlaying Image Axes\n", - "\n", - "For the common use case of plotting image axis on an `ArrayDisplay`, the `set_line_hillas()` method is provided for convenience. The following example shows its use: " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "cb0b242c-4fc4-48ff-b8f2-2d57b852deea", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "from astropy.coordinates import SkyCoord\n", - "from ctapipe.calib import CameraCalibrator\n", - "from ctapipe.image import ImageProcessor\n", - "from ctapipe.io import EventSource\n", - "from ctapipe.reco import ShowerProcessor\n", - "from ctapipe.utils import get_dataset_path\n", - "from ctapipe.visualization import ArrayDisplay\n", - "from IPython import display\n", - "from matplotlib.animation import FuncAnimation\n", - "\n", - "input_url = \"dataset://gamma_LaPalma_baseline_20Zd_180Az_prod3b_test.simtel.gz\"" - ] - }, - { - "cell_type": "markdown", - "id": "f073a3a5-5db3-43fa-8734-deea0de7b7b4", - "metadata": { - "tags": [] - }, - "source": [ - "First, we define a function to plot the array with overlaid lines for the image axes" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b24b2c6b-60fb-4a5c-b433-9be76ab26cee", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "def plot_event(event, subarray, ax):\n", - " \"\"\"\n", - " Draw an ArrayDisplay with image axes and the\n", - " true and reconstructed impact position overlaid\n", - " \"\"\"\n", - "\n", - " array_pointing = SkyCoord(\n", - " az=event.pointing.array_azimuth,\n", - " alt=event.pointing.array_altitude,\n", - " frame=\"altaz\",\n", - " )\n", - "\n", - " angle_offset = event.pointing.array_azimuth\n", - " disp = ArrayDisplay(subarray, axes=ax)\n", - "\n", - " hillas_dict = {tid: tel.parameters.hillas for tid, tel in event.dl1.tel.items()}\n", - " core_dict = {tid: tel.parameters.core.psi for tid, tel in event.dl1.tel.items()}\n", - "\n", - " disp.set_line_hillas(\n", - " hillas_dict,\n", - " core_dict,\n", - " 500,\n", - " )\n", - "\n", - " reco_shower = event.dl2.stereo.geometry[\"HillasReconstructor\"]\n", - "\n", - " ax.scatter(\n", - " event.simulation.shower.core_x,\n", - " event.simulation.shower.core_y,\n", - " s=200,\n", - " c=\"k\",\n", - " marker=\"x\",\n", - " label=\"True Impact\",\n", - " )\n", - " ax.scatter(\n", - " reco_shower.core_x,\n", - " reco_shower.core_y,\n", - " s=200,\n", - " c=\"r\",\n", - " marker=\"x\",\n", - " label=\"Estimated Impact\",\n", - " )\n", - "\n", - " ax.legend()" - ] - }, - { - "cell_type": "markdown", - "id": "47c2e70a-f1f3-4922-92e6-12fffe0566a4", - "metadata": {}, - "source": [ - "Now, we can loop through some events and plot them. Here we apply default calibration, image processing, and reconstruction, however it is better to use `ctapipe-process` with a well-defined configuration to do this in reality. Note that some events will not have images bright enough to do parameterization or reconstruction, so they will have no image axis lines or no estimated impact position." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "502fe577-3e8a-4e6d-90f3-8db8fbb57b90", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(5, 3, figsize=(20, 40), constrained_layout=True)\n", - "ax = ax.ravel()\n", - "\n", - "with EventSource(input_url, max_events=15, focal_length_choice=\"EQUIVALENT\") as source:\n", - " calib = CameraCalibrator(subarray=source.subarray)\n", - " process_images = ImageProcessor(subarray=source.subarray)\n", - " process_shower = ShowerProcessor(subarray=source.subarray)\n", - "\n", - " for i, event in enumerate(source):\n", - " calib(event)\n", - " process_images(event)\n", - " process_shower(event)\n", - " plot_event(event, source.subarray, ax=ax[i])" - ] - } - ], - "metadata": { - "language_info": { - "codemirror_mode": { - "name": "ipython" - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/docs/user-guide/examples/camera_display.ipynb b/docs/user-guide/examples/camera_display.ipynb deleted file mode 100644 index 3d1274ee515..00000000000 --- a/docs/user-guide/examples/camera_display.ipynb +++ /dev/null @@ -1,560 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Displaying Camera Images" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import astropy.coordinates as c\n", - "import astropy.units as u\n", - "import matplotlib.pylab as plt\n", - "import numpy as np\n", - "from ctapipe.coordinates import CameraFrame, EngineeringCameraFrame, TelescopeFrame\n", - "from ctapipe.image import hillas_parameters, tailcuts_clean, toymodel\n", - "from ctapipe.instrument import SubarrayDescription\n", - "from ctapipe.visualization import CameraDisplay" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "First, let's create a fake Cherenkov image from a given `CameraGeometry` and fill it with some data that we can draw later." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# load an example camera geometry from a simulation file\n", - "subarray = SubarrayDescription.read(\"dataset://gamma_prod5.simtel.zst\")\n", - "geom = subarray.tel[100].camera.geometry\n", - "\n", - "# create a fake camera image to display:\n", - "model = toymodel.Gaussian(\n", - " x=0.2 * u.m,\n", - " y=0.0 * u.m,\n", - " width=0.05 * u.m,\n", - " length=0.15 * u.m,\n", - " psi=\"35d\",\n", - ")\n", - "\n", - "image, sig, bg = model.generate_image(geom, intensity=1500, nsb_level_pe=10)\n", - "mask = tailcuts_clean(geom, image, picture_thresh=15, boundary_thresh=5)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "geom" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Displaying Images\n", - "\n", - "The simplest plot is just to generate a CameraDisplay with an image in its constructor. A figure and axis will be created automatically" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "CameraDisplay(geom)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "You can also specify the initial `image`, `cmap` and `norm` (colomap and normalization, see below), `title` to use. You can specify `ax` if you want to draw the camera on an existing *matplotlib* `Axes` object (otherwise one is created).\n", - "\n", - "To change other options, or to change options dynamically, you can call the relevant functions of the `CameraDisplay` object that is returned. For example to add a color bar, call `add_colorbar()`, or to change the color scale, modify the `cmap` or `norm` properties directly. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Choosing a coordinate frame\n", - "\n", - "The `CameraGeometry` object contains a `ctapipe.coordinates.Frame` used by `CameraDisplay` to draw the camera in the correct orientation and distance units. The default frame is the `CameraFrame`, which will display the camera in units of *meters* and with an orientation that the top of the camera (when parked) is aligned to the X-axis. To show the camera in another orientation, it's useful to apply a coordinate transform to the `CameraGeometry` before passing it to the `CameraDisplay`. The following `Frames` are supported:\n", - "* `EngineeringCameraFrame` : similar to CameraFrame, but with the top of the camera aligned to the Y axis\n", - "* `TelescopeFrame`: In *degrees* (on the sky) coordinates relative to the telescope Alt/Az pointing position, with the Alt axis pointing upward. \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(1, 3, figsize=(15, 4))\n", - "CameraDisplay(geom, image=image, ax=ax[0])\n", - "CameraDisplay(geom.transform_to(EngineeringCameraFrame()), image=image, ax=ax[1])\n", - "CameraDisplay(geom.transform_to(TelescopeFrame()), image=image, ax=ax[2])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note the the name of the Frame appears in the lower-right corner" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "For the rest of this demo, let's use the `TelescopeFrame`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "geom_camframe = geom\n", - "geom = geom_camframe.transform_to(EngineeringCameraFrame())" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Changing the color map and scale" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "CameraDisplay supports any [matplotlib color map](https://matplotlib.org/stable/tutorials/colors/colormaps.html)\n", - "It is **highly recommended** to use a *perceptually uniform* map, unless you have a good reason not to." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(1, 3, figsize=(15, 4))\n", - "for ii, cmap in enumerate([\"PuOr_r\", \"rainbow\", \"twilight\"]):\n", - " disp = CameraDisplay(geom, image=image, ax=ax[ii], title=cmap)\n", - " disp.add_colorbar()\n", - " disp.cmap = cmap" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "By default the minimum and maximum of the color bar are set automatically by the data in the image. To choose fixed limits, use:`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(1, 3, figsize=(15, 4))\n", - "for ii, minmax in enumerate([(10, 50), (-10, 10), (1, 100)]):\n", - " disp = CameraDisplay(geom, image=image, ax=ax[ii], title=minmax)\n", - " disp.add_colorbar()\n", - " disp.set_limits_minmax(minmax[0], minmax[1])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Or you can set the maximum limit by percentile of the charge distribution:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(1, 3, figsize=(15, 4))\n", - "for ii, pct in enumerate([30, 50, 90]):\n", - " disp = CameraDisplay(geom, image=image, ax=ax[ii], title=f\"{pct} %\")\n", - " disp.add_colorbar()\n", - " disp.set_limits_percent(pct)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Using different normalizations" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "You can choose from several preset normalizations (lin, log, symlog) and also provide a custom normalization, for example a `PowerNorm`:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from matplotlib.colors import PowerNorm\n", - "\n", - "fig, axes = plt.subplots(2, 2, figsize=(14, 10))\n", - "norms = [\"lin\", \"log\", \"symlog\", PowerNorm(0.5)]\n", - "\n", - "for norm, ax in zip(norms, axes.flatten()):\n", - " disp = CameraDisplay(geom, image=image, ax=ax)\n", - " disp.norm = norm\n", - " disp.add_colorbar()\n", - " ax.set_title(str(norm))\n", - "\n", - "axes[1, 1].set_title(\"PowerNorm(0.5)\")\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Overlays" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Marking pixels\n", - "\n", - "here we will mark pixels in the image mask. That will change their outline color" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(1, 2, figsize=(10, 4))\n", - "disp = CameraDisplay(\n", - " geom, image=image, cmap=\"gray\", ax=ax[0], title=\"Image mask in green\"\n", - ")\n", - "disp.highlight_pixels(mask, alpha=0.8, linewidth=2, color=\"green\")\n", - "\n", - "disp = CameraDisplay(\n", - " geom, image=image, cmap=\"gray\", ax=ax[1], title=\"Image mask in green (zoom)\"\n", - ")\n", - "disp.highlight_pixels(mask, alpha=1, linewidth=3, color=\"green\")\n", - "\n", - "ax[1].set_ylim(-0.5, 0.5)\n", - "ax[1].set_xlim(-0.5, 0.5)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Drawing a Hillas-parameter ellipse" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "For this, we will first compute some Hillas Parameters in the current frame:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "clean_image = image.copy()\n", - "clean_image[~mask] = 0\n", - "hillas = hillas_parameters(geom, clean_image)\n", - "\n", - "plt.figure(figsize=(6, 6))\n", - "disp = CameraDisplay(geom, image=image, cmap=\"gray_r\")\n", - "disp.highlight_pixels(mask, alpha=0.5, color=\"dodgerblue\")\n", - "disp.overlay_moments(hillas, color=\"red\", linewidth=3, with_label=False)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Drawing a marker at a coordinate\n", - "\n", - "This depends on the coordinate frame of the `CameraGeometry`. Here we will sepcify the coordinate the `EngineerngCameraFrame`, but if you have enough information to do the coordinate transform, you could use `ICRS` coordinates and overlay star positions. `CameraDisplay` will convert the coordinate you pass in to the `Frame` of the display automatically (if sufficient frame attributes are set). \n", - "\n", - "Note that the parameter `keep_old` is False by default, meaning adding a new point will clear the previous ones (useful for animations, but perhaps unexpected for a static plot). Set it to `True` to plot multiple markers." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.figure(figsize=(6, 6))\n", - "disp = CameraDisplay(geom, image=image, cmap=\"gray_r\")\n", - "\n", - "coord = c.SkyCoord(x=0.5 * u.m, y=0.7 * u.m, frame=geom.frame)\n", - "coord_in_another_frame = c.SkyCoord(x=0.5 * u.m, y=0.7 * u.m, frame=CameraFrame())\n", - "disp.overlay_coordinate(coord, markersize=20, marker=\"*\")\n", - "disp.overlay_coordinate(\n", - " coord_in_another_frame, markersize=20, marker=\"*\", keep_old=True\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Generating an animation\n", - "\n", - "Here we will make an animation of fake events by re-using a single display (much faster than generating a new one each time) " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from IPython import display\n", - "from matplotlib.animation import FuncAnimation\n", - "\n", - "subarray = SubarrayDescription.read(\"dataset://gamma_prod5.simtel.zst\")\n", - "geom = subarray.tel[1].camera.geometry\n", - "\n", - "fov = 1.0\n", - "maxwid = 0.05\n", - "maxlen = 0.1\n", - "\n", - "fig, ax = plt.subplots(1, 1, figsize=(8, 6))\n", - "disp = CameraDisplay(geom, ax=ax) # we only need one display (it can be re-used)\n", - "disp.cmap = \"inferno\"\n", - "disp.add_colorbar(ax=ax)\n", - "\n", - "\n", - "def update(frame):\n", - " \"\"\"this function will be called for each frame of the animation\"\"\"\n", - " x, y = np.random.uniform(-fov, fov, size=2)\n", - " width = np.random.uniform(0.01, maxwid)\n", - " length = np.random.uniform(width, maxlen)\n", - " angle = np.random.uniform(0, 180)\n", - " intens = width * length * (5e4 + 1e5 * np.random.exponential(2))\n", - "\n", - " model = toymodel.Gaussian(\n", - " x=x * u.m,\n", - " y=y * u.m,\n", - " width=width * u.m,\n", - " length=length * u.m,\n", - " psi=angle * u.deg,\n", - " )\n", - " image, _, _ = model.generate_image(\n", - " geom,\n", - " intensity=intens,\n", - " nsb_level_pe=5,\n", - " )\n", - " disp.image = image\n", - "\n", - "\n", - "# Create the animation and convert to a displayable video:\n", - "anim = FuncAnimation(fig, func=update, frames=10, interval=200)\n", - "plt.close(fig) # so it doesn't display here\n", - "video = anim.to_html5_video()\n", - "display.display(display.HTML(video))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Using CameraDisplays interactively\n", - "\n", - "`CameraDisplays` can be used interactivly whe displayed in a window, and also when using Jupyter notebooks/lab with appropriate backends. \n", - "\n", - "When this is the case, the same `CameraDisplay` object can be re-used. We can't show this here in the documentation, but creating an animation when in a matplotlib window is quite easy! Try this in an interactive ipython session:\n", - "\n", - "### Running interactive displays in a matplotlib window\n", - "\n", - "```sh\n", - "ipython -i --maplotlib=auto\n", - "```\n", - "\n", - "That will open an ipython session with matplotlib graphics in a separate thread, meaning that you can type code and interact with plots simultaneneously. \n", - "\n", - "In the ipython session try running the following code and you will see an animation (here in the documentation, it will of course be static)\n", - "\n", - "First we load some real data so we have a nice image to view:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "from ctapipe.io import EventSource\n", - "from ctapipe.visualization import CameraDisplay\n", - "import numpy as np\n", - "\n", - "DATA = \"dataset://gamma_20deg_0deg_run1___cta-prod5-lapalma_desert-2158m-LaPalma-dark_100evts.simtel.zst\"\n", - "\n", - "with EventSource(\n", - " DATA,\n", - " max_events=1,\n", - " focal_length_choice=\"EQUIVALENT\",\n", - ") as source:\n", - " event = next(iter(source))\n", - "\n", - "tel_id = list(event.r0.tel.keys())[0]\n", - "geom = source.subarray.tel[tel_id].camera.geometry\n", - "waveform = event.r0.tel[tel_id].waveform\n", - "n_chan, n_pix, n_samp = waveform.shape" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Running the following the will bring up a window and animate the shower image as a function of time. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "disp = CameraDisplay(geom)\n", - "\n", - "for ii in range(n_samp):\n", - " disp.image = waveform[0, :, ii]\n", - " plt.pause(0.1) # this lets matplotlib re-draw the scene" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The output will be similar to the static animation created as follows:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(1, 1)\n", - "disp = CameraDisplay(geom, ax=ax)\n", - "disp.add_colorbar()\n", - "disp.autoscale = False\n", - "\n", - "\n", - "def draw_sample(frame):\n", - " ax.set_title(f\"sample: {frame}\")\n", - " disp.set_limits_minmax(200, 400)\n", - " disp.image = waveform[0, :, frame]\n", - "\n", - "\n", - "anim = FuncAnimation(fig, func=draw_sample, frames=n_samp, interval=100)\n", - "plt.close(fig) # so it doesn't display here\n", - "video = anim.to_html5_video()\n", - "display.display(display.HTML(video))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Making it clickable \n", - "\n", - "Also when running in a window, you can enable the `disp.enable_pixel_picker()` option. This will then allow the user to click a pixel and a function will run. By default the function simply prints the pixel and value to stdout, however you can override the function `on_pixel_clicked(pix_id)` to do anything you want by making a subclass" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "class MyCameraDisplay(CameraDisplay):\n", - " def on_pixel_clicked(self, pix_id):\n", - " print(f\"{pix_id=} has value {self.image[pix_id]:.2f}\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "disp = MyCameraDisplay(geom, image=image)\n", - "disp.enable_pixel_picker()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "then, when a user clicks a pixel it would print:\n", - "```\n", - "pixel 5 has value 2.44\n", - "```" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "anaconda-cloud": {}, - "language_info": { - "codemirror_mode": { - "name": "ipython" - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/docs/user-guide/examples/containers.ipynb b/docs/user-guide/examples/containers.ipynb deleted file mode 100644 index 603139f6f41..00000000000 --- a/docs/user-guide/examples/containers.ipynb +++ /dev/null @@ -1,396 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Using Container classes\n", - "\n", - "`ctapipe.core.Container` is the base class for all event-wise data classes in ctapipe. It works like a object-relational mapper, in that it defines a set of `Fields` along with their metadata (description, unit, default), which can be later translated automatially into an output table using a `ctapipe.io.TableWriter`." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.core import Container, Field, Map\n", - "import numpy as np\n", - "from astropy import units as u\n", - "from functools import partial" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's define a few example containers with some dummy fields in them:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "class SubContainer(Container):\n", - " junk = Field(-1, \"Some junk\")\n", - " value = Field(0.0, \"some value\", unit=u.deg)\n", - "\n", - "\n", - "class TelContainer(Container):\n", - " # defaults should match the other requirements, e.g. the defaults\n", - " # should have the correct unit. It most often also makes sense to use\n", - " # an invalid value marker like nan for floats or -1 for positive integers\n", - " # as default\n", - " tel_id = Field(-1, \"telescope ID number\")\n", - "\n", - " # For mutable structures like lists, arrays or containers, use a `default_factory` function or class\n", - " # not an instance to assure each container gets a fresh instance and there is no hidden\n", - " # shared state between containers.\n", - " image = Field(default_factory=lambda: np.zeros(10), description=\"camera pixel data\")\n", - "\n", - "\n", - "class EventContainer(Container):\n", - " event_id = Field(-1, \"event id number\")\n", - "\n", - " tels_with_data = Field(\n", - " default_factory=list, description=\"list of telescopes with data\"\n", - " )\n", - " sub = Field(\n", - " default_factory=SubContainer, description=\"stuff\"\n", - " ) # a sub-container in the hierarchy\n", - "\n", - " # A Map is like a defaultdictionary with a specific container type as default.\n", - " # This can be used to e.g. store a container per telescope\n", - " # we use partial here to automatically get a function that creates a map with the correct container type\n", - " # as default\n", - " tel = Field(default_factory=partial(Map, TelContainer), description=\"telescopes\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Basic features" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ev = EventContainer()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Check that default values are automatically filled in" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(ev.event_id)\n", - "print(ev.sub)\n", - "print(ev.tel)\n", - "print(ev.tel.keys())\n", - "\n", - "# default dict access will create container:\n", - "print(ev.tel[1])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "print the dict representation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(ev)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We also get docstrings \"for free\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "?EventContainer" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "?SubContainer" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "values can be set as normal for a class:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ev.event_id = 100\n", - "ev.event_id" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ev.as_dict() # by default only shows the bare items, not sub-containers (See later)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ev.as_dict(recursive=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "and we can add a few of these to the parent container inside the tel dict:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ev.tel[10] = TelContainer()\n", - "ev.tel[5] = TelContainer()\n", - "ev.tel[42] = TelContainer()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# because we are using a default_factory to handle mutable defaults, the images are actually different:\n", - "ev.tel[42].image is ev.tel[32]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Be careful to use the `default_factory` mechanism for mutable fields, see this **negative** example:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "class DangerousContainer(Container):\n", - " image = Field(\n", - " np.zeros(10),\n", - " description=\"Attention!!!! Globally mutable shared state. Use default_factory instead\",\n", - " )\n", - "\n", - "\n", - "c1 = DangerousContainer()\n", - "c2 = DangerousContainer()\n", - "\n", - "c1.image[5] = 9999\n", - "\n", - "print(c1.image)\n", - "print(c2.image)\n", - "print(c1.image is c2.image)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ev.tel" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Converion to dictionaries" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ev.as_dict()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ev.as_dict(recursive=True, flatten=False)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "for serialization to a table, we can even flatten the output into a single set of columns" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ev.as_dict(recursive=True, flatten=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setting and clearing values" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ev.tel[5].image[:] = 9\n", - "print(ev)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ev.reset()\n", - "ev.as_dict(recursive=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## look at a pre-defined Container" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.containers import SimulatedShowerContainer" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "?SimulatedShowerContainer" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "shower = SimulatedShowerContainer()\n", - "shower" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Container prefixes\n", - "\n", - "To store the same container in the same table in a file or give more information, containers support setting\n", - "a custom prefix:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "c1 = SubContainer(junk=5, value=3, prefix=\"foo\")\n", - "c2 = SubContainer(junk=10, value=9001, prefix=\"bar\")\n", - "\n", - "# create a common dict with data from both containers:\n", - "d = c1.as_dict(add_prefix=True)\n", - "d.update(c2.as_dict(add_prefix=True))\n", - "d" - ] - } - ], - "metadata": { - "anaconda-cloud": {}, - "language_info": { - "codemirror_mode": { - "name": "ipython" - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/docs/user-guide/examples/convert_images_to_2d.ipynb b/docs/user-guide/examples/convert_images_to_2d.ipynb deleted file mode 100644 index a38e135ab2d..00000000000 --- a/docs/user-guide/examples/convert_images_to_2d.ipynb +++ /dev/null @@ -1,208 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Convert camera images to pixels on a s square grid" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.utils import get_dataset_path\n", - "from ctapipe.visualization import CameraDisplay\n", - "from ctapipe.instrument import SubarrayDescription\n", - "from ctapipe.io import EventSource\n", - "from ctapipe.image.toymodel import Gaussian\n", - "import matplotlib.pyplot as plt\n", - "import astropy.units as u" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# get the subarray from an example file\n", - "subarray = SubarrayDescription.read(\"dataset://gamma_prod5.simtel.zst\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Geometries with square pixels\n", - "\n", - "Define a camera geometry and generate a dummy image:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "geom = subarray.tel[40].camera.geometry\n", - "model = Gaussian(\n", - " x=0.05 * u.m,\n", - " y=0.05 * u.m,\n", - " width=0.01 * u.m,\n", - " length=0.05 * u.m,\n", - " psi=\"30d\",\n", - ")\n", - "_, image, _ = model.generate_image(geom, intensity=500, nsb_level_pe=3)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "CameraDisplay(geom, image)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The `CameraGeometry` has functions to convert the 1d image arrays to 2d arrays and back to the 1d array:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "image_square = geom.image_to_cartesian_representation(image)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.imshow(image_square)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "image_1d = geom.image_from_cartesian_representation(image_square)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "CameraDisplay(geom, image_1d)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Geometries with hexagonal pixels\n", - "\n", - "Define a camera geometry and generate a dummy image:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "geom = subarray.tel[1].camera.geometry\n", - "model = Gaussian(\n", - " x=0.5 * u.m,\n", - " y=0.5 * u.m,\n", - " width=0.1 * u.m,\n", - " length=0.2 * u.m,\n", - " psi=\"30d\",\n", - ")\n", - "_, image, _ = model.generate_image(geom, intensity=5000)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "CameraDisplay(geom, image)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "image_square = geom.image_to_cartesian_representation(image)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Conversion into square geometry\n", - "\n", - "Since the resulting array has square pixels, the pixel grid has to be rotated and distorted.\n", - "This is reversible (The `image_from_cartesian_representation` method takes care of this):" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.imshow(image_square)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "image_1d = geom.image_from_cartesian_representation(image_square)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "disp = CameraDisplay(geom, image_1d)" - ] - } - ], - "metadata": { - "language_info": { - "codemirror_mode": { - "name": "ipython" - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/docs/user-guide/examples/dilate_image.ipynb b/docs/user-guide/examples/dilate_image.ipynb deleted file mode 100644 index e833a07f142..00000000000 --- a/docs/user-guide/examples/dilate_image.ipynb +++ /dev/null @@ -1,128 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Basic Image Cleaning and Dilation\n", - "\n", - "Here we create an example shower image, do a tail-cuts (picture/boundary) cleaning, and then dilate the resulting cleaning mask by several neighbor pixels" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib inline\n", - "from matplotlib import pyplot as plt\n", - "import astropy.units as u" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.instrument import SubarrayDescription\n", - "from ctapipe.visualization import CameraDisplay\n", - "from ctapipe.image import toymodel, tailcuts_clean, dilate\n", - "\n", - "# Load a camera from an example file\n", - "subarray = SubarrayDescription.read(\"dataset://gamma_prod5.simtel.zst\")\n", - "geom = subarray.tel[100].camera.geometry\n", - "\n", - "# Create a fake camera image to display:\n", - "model = toymodel.Gaussian(\n", - " x=0.2 * u.m,\n", - " y=0.0 * u.m,\n", - " width=0.05 * u.m,\n", - " length=0.15 * u.m,\n", - " psi=\"35d\",\n", - ")\n", - "\n", - "image, sig, bg = model.generate_image(geom, intensity=1500, nsb_level_pe=5)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Apply the image cleaning:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cleanmask = tailcuts_clean(geom, image, picture_thresh=10, boundary_thresh=5)\n", - "clean = image.copy()\n", - "clean[~cleanmask] = 0.0\n", - "\n", - "disp = CameraDisplay(geom, image=image)\n", - "disp.highlight_pixels(cleanmask, color=\"red\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now dialte the mask a few times:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.image.cleaning import dilate" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def show_dilate(mask, times=1):\n", - " m = mask.copy()\n", - " for ii in range(times):\n", - " m = dilate(geom, m)\n", - " CameraDisplay(\n", - " geom, image=(m.astype(int) + mask.astype(int)), title=\"dilate{}\".format(times)\n", - " )" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.figure(figsize=(18, 3))\n", - "\n", - "for ii in range(0, 6):\n", - " plt.subplot(1, 7, ii + 1)\n", - " show_dilate(cleanmask.copy(), times=ii)" - ] - } - ], - "metadata": { - "language_info": { - "codemirror_mode": { - "name": "ipython" - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/docs/user-guide/examples/index.rst b/docs/user-guide/examples/index.rst deleted file mode 100644 index 71eb0c33078..00000000000 --- a/docs/user-guide/examples/index.rst +++ /dev/null @@ -1,36 +0,0 @@ -.. _examples: - -======== -Examples -======== - -Some lower-level examples of features of the ctapipe API (see -the Tutorials section for more complete examples) - - -.. toctree:: - :maxdepth: 1 - :caption: Visualization - - camera_display - array_display - - -.. toctree:: - :maxdepth: 1 - :caption: Algorithms - - dilate_image - nd_interpolation - convert_images_to_2d - - -.. toctree:: - :maxdepth: 1 - :caption: Core functionality - - InstrumentDescription - containers - Tools - provenance - table_writer_reader diff --git a/docs/user-guide/examples/nd_interpolation.ipynb b/docs/user-guide/examples/nd_interpolation.ipynb deleted file mode 100644 index f98d71c8cef..00000000000 --- a/docs/user-guide/examples/nd_interpolation.ipynb +++ /dev/null @@ -1,207 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Use N-dimensional Histogram functionality and Interpolation\n", - " \n", - "* could be used for example to read and interpolate an lookup table or IRF.\n", - "* In this example, we load a sample energy reconstruction lookup-table from a FITS file\n", - "* In this case it is only in 2D cube (to keep the file size small): `SIZE` vs `IMPACT-DISTANCE`, however the same method will work for any dimensionality" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pylab as plt\n", - "import numpy as np\n", - "from astropy.io import fits\n", - "from scipy.interpolate import RegularGridInterpolator\n", - "\n", - "from ctapipe.utils.datasets import get_dataset_path\n", - "from ctapipe.utils import Histogram\n", - "\n", - "%matplotlib inline" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### load an example datacube \n", - "(an energy table generated for a\n", - "small subset of HESS simulations) to use as a lookup table. Here\n", - "we will use the `Histogram` class, which automatically loads both\n", - "the data cube and creates arrays for the coordinates of each\n", - "axis." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "testfile = get_dataset_path(\"hess_ct001_energylut.fits.gz\")\n", - "energy_hdu = fits.open(testfile)[\"MEAN\"]\n", - "energy_table = Histogram.from_fits(energy_hdu)\n", - "print(energy_table)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### construct an interpolator that we can use to get values at any point:\n", - "\n", - "Here we will use a `RegularGridInterpolator`, since it is the most appropriate for this type of data, but others are available (see the SciPy documentation)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "centers = [energy_table.bin_centers(ii) for ii in range(energy_table.ndims)]\n", - "energy_lookup = RegularGridInterpolator(\n", - " centers, energy_table.hist, bounds_error=False, fill_value=-100\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "`energy_lookup` is now just a continuous function of `log(SIZE)`,\n", - "`DISTANCE` in m. \n", - "\n", - "### Now plot some curves from the interpolator. \n", - "\n", - "Note that the LUT we used is does not have very high statistics,\n", - "so the interpolation starts to be affected by noise at the high\n", - "end. In a real case, we would want to use a table that has been\n", - "sanitized (smoothed and extrapolated)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "lsize = np.linspace(1.5, 5.0, 100)\n", - "dists = np.linspace(50, 100, 5)\n", - "\n", - "plt.figure()\n", - "plt.title(\"Variation of energy with size and impact distance\")\n", - "plt.xlabel(\"SIZE (P.E.)\")\n", - "plt.ylabel(\"ENERGY (TeV)\")\n", - "\n", - "for dist in dists:\n", - " plt.plot(\n", - " 10**lsize,\n", - " 10 ** energy_lookup((lsize, dist)),\n", - " \"+-\",\n", - " label=\"DIST={:.1f} m\".format(dist),\n", - " )\n", - "\n", - "plt.legend(loc=\"best\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Using the interpolator, reinterpolate the lookup table onto an $N \\times N$\n", - "grid (regardless of its original dimensions):" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "N = 300\n", - "xmin, xmax = energy_table.bin_centers(0)[0], energy_table.bin_centers(0)[-1]\n", - "ymin, ymax = energy_table.bin_centers(1)[0], energy_table.bin_centers(1)[-1]\n", - "xx, yy = np.linspace(xmin, xmax, N), np.linspace(ymin, ymax, N)\n", - "X, Y = np.meshgrid(xx, yy)\n", - "points = list(zip(X.ravel(), Y.ravel()))\n", - "E = energy_lookup(points).reshape((N, N))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now, let's plot the original table and the new one (E). The color bar shows $\\log_{10}(E)$ in TeV" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.figure(figsize=(12, 5))\n", - "plt.nipy_spectral()\n", - "\n", - "# the uninterpolated table\n", - "plt.subplot(1, 2, 1)\n", - "plt.xlim(1.5, 5)\n", - "plt.ylim(0, 500)\n", - "plt.xlabel(\"log10(SIZE)\")\n", - "plt.ylabel(\"Impact Dist (m)\")\n", - "plt.pcolormesh(\n", - " energy_table.bin_centers(0), energy_table.bin_centers(1), energy_table.hist.T\n", - ")\n", - "plt.title(\"Raw table, uninterpolated {0}\".format(energy_table.hist.T.shape))\n", - "cb = plt.colorbar()\n", - "cb.set_label(\"$\\log_{10}(E/\\mathrm{TeV})$\")\n", - "\n", - "# the interpolated table\n", - "plt.subplot(1, 2, 2)\n", - "plt.pcolormesh(np.linspace(xmin, xmax, N), np.linspace(ymin, ymax, N), E)\n", - "plt.xlim(1.5, 5)\n", - "plt.ylim(0, 500)\n", - "plt.xlabel(\"log10(SIZE)\")\n", - "plt.ylabel(\"Impact Dist(m)\")\n", - "plt.title(\"Interpolated to a ({0}, {0}) grid\".format(N))\n", - "cb = plt.colorbar()\n", - "cb.set_label(\"$\\log_{10}(E/\\mathrm{TeV})$\")\n", - "\n", - "plt.tight_layout()\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In the high-stats central region, we get a nice smooth interpolation function. Of course we can see that there are a few more steps to take before using this table:\n", - "* need to deal with cases where the table had low stats near the edges (smooth or extrapolate, or set bounds)\n", - "* may need to smooth the table even where there are sufficient stats, to avoid wiggles" - ] - } - ], - "metadata": { - "anaconda-cloud": {}, - "language_info": { - "codemirror_mode": { - "name": "ipython" - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/docs/user-guide/examples/provenance.ipynb b/docs/user-guide/examples/provenance.ipynb deleted file mode 100644 index 2743eed9261..00000000000 --- a/docs/user-guide/examples/provenance.ipynb +++ /dev/null @@ -1,199 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Using the ctapipe Provenance service\n", - "\n", - "The provenance functionality is used automatically when you use most of ctapipe functionality (particularly `ctapipe.core.Tool` and functions in `ctapipe.io` and `ctapipe.utils`), so normally you don't have to work with it directly. It tracks both input and output files, as well as details of the machine and software environment on which a Tool executed. \n", - "\n", - "Here we show some very low-level functions of this system:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.core import Provenance\n", - "from pprint import pprint" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Activities\n", - "\n", - "The basis of Provenance is an *activity*, which is generally an executable or step in a script. Activities can be nested (e.g. with sub-activities), as shown below, but normally this is not required:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "p = Provenance() # note this is a singleton, so only ever one global provenence object\n", - "p.clear()\n", - "p.start_activity()\n", - "p.add_input_file(\"test.txt\")\n", - "\n", - "p.start_activity(\"sub\")\n", - "p.add_input_file(\"subinput.txt\")\n", - "p.add_input_file(\"anothersubinput.txt\")\n", - "p.add_output_file(\"suboutput.txt\")\n", - "p.finish_activity(\"sub\")\n", - "\n", - "p.start_activity(\"sub2\")\n", - "p.add_input_file(\"sub2input.txt\")\n", - "p.finish_activity(\"sub2\")\n", - "\n", - "p.finish_activity()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "p.finished_activity_names" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Activities have associated input and output *entities* (files or other objects)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "[(x[\"activity_name\"], x[\"input\"]) for x in p.provenance]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Activities track when they were started and finished:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "[(x[\"activity_name\"], x[\"duration_min\"]) for x in p.provenance]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Full provenance\n", - "\n", - "The provence object is a list of activitites, and for each lots of details are collected:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "p.provenance[0]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This can be better represented in JSON:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(p.as_json(indent=2))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Storing provenance info in output files\n", - "\n", - "* already this can be stored in something like an HDF5 file header, which allows hierarchies.\n", - "* Try to flatted the data so it can be stored in a key=value header in a **FITS file** (using the FITS extended keyword convention to allow >8 character keywords), or as a table " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def flatten_dict(y):\n", - " out = {}\n", - "\n", - " def flatten(x, name=\"\"):\n", - " if type(x) is dict:\n", - " for a in x:\n", - " flatten(x[a], name + a + \".\")\n", - " elif type(x) is list:\n", - " i = 0\n", - " for a in x:\n", - " flatten(a, name + str(i) + \".\")\n", - " i += 1\n", - " else:\n", - " out[name[:-1]] = x\n", - "\n", - " flatten(y)\n", - " return out" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "d = dict(activity=p.provenance)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "pprint(flatten_dict(d))" - ] - } - ], - "metadata": { - "language_info": { - "codemirror_mode": { - "name": "ipython" - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/docs/user-guide/examples/table_writer_reader.ipynb b/docs/user-guide/examples/table_writer_reader.ipynb deleted file mode 100644 index aa27756d0c4..00000000000 --- a/docs/user-guide/examples/table_writer_reader.ipynb +++ /dev/null @@ -1,456 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Writing Containers to a tabular format\n", - "\n", - "The `TableWriter`/`TableReader` sub-classes allow you to write a `ctapipe.core.Container` class and its meta-data to an output table. They treat the `Field`s in the `Container` as columns in the output, and automatically generate a schema. Here we will go through an example of writing out data and reading it back with *Pandas*, *PyTables*, and a `ctapipe.io.TableReader`:\n", - "\n", - "In this example, we will use the `HDF5TableWriter`, which writes to HDF5 datasets using *PyTables*. Currently this is the only implemented TableWriter." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Caveats to think about:\n", - "* vector columns in Containers *can* be written, but some lilbraries like Pandas can not read those (so you must use pytables or astropy to read outputs that have vector columns)\n", - "* units are stored in the table metadata, but some libraries like Pandas ignore them and all other metadata" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Create some example Containers" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.io import HDF5TableWriter\n", - "from ctapipe.core import Container, Field\n", - "from astropy import units as u\n", - "import numpy as np" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "class VariousTypesContainer(Container):\n", - "\n", - " a_int = Field(int, \"some int value\")\n", - " a_float = Field(float, \"some float value with a unit\", unit=u.m)\n", - " a_bool = Field(bool, \"some bool value\")\n", - " a_np_int = Field(np.int64, \"a numpy int\")\n", - " a_np_float = Field(np.float64, \"a numpy float\")\n", - " a_np_bool = Field(np.bool_, \"np.bool\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "let's also make a dummy stream (generator) that will create a series of these containers" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def create_stream(n_event):\n", - "\n", - " data = VariousTypesContainer()\n", - " for i in range(n_event):\n", - "\n", - " data.a_int = int(i)\n", - " data.a_float = float(i) * u.cm # note unit conversion will happen\n", - " data.a_bool = (i % 2) == 0\n", - " data.a_np_int = np.int64(i)\n", - " data.a_np_float = np.float64(i)\n", - " data.a_np_bool = np.bool_((i % 2) == 0)\n", - "\n", - " yield data" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for data in create_stream(2):\n", - "\n", - " for key, val in data.items():\n", - "\n", - " print(\"{}: {}, type : {}\".format(key, val, type(val)))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Writing the Data (and good practices)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Always use context managers with IO classes, as they will make sure the underlying resources are properly closed in case of errors:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "try:\n", - " with HDF5TableWriter(\"container.h5\", group_name=\"data\") as h5_table:\n", - "\n", - " for data in create_stream(10):\n", - "\n", - " h5_table.write(\"table\", data)\n", - " 0 / 0\n", - "except Exception as err:\n", - " print(\"FAILED:\", err)\n", - "print(\"Done\")\n", - "\n", - "h5_table.h5file.isopen == False" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "!ls container.h5" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Appending new Containers" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To append some new containers we need to set the writing in append mode by using: 'mode=a'. But let's now first look at what happens if we don't." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for i in range(2):\n", - "\n", - " with HDF5TableWriter(\n", - " \"container.h5\", mode=\"w\", group_name=\"data_{}\".format(i)\n", - " ) as h5_table:\n", - "\n", - " for data in create_stream(10):\n", - "\n", - " h5_table.write(\"table\", data)\n", - "\n", - " print(h5_table.h5file)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "!rm -f container.h5" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Ok so the writer destroyed the content of the file each time it opens the file. Now let's try to append some data group to it! (using mode='a')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for i in range(2):\n", - "\n", - " with HDF5TableWriter(\n", - " \"container.h5\", mode=\"a\", group_name=\"data_{}\".format(i)\n", - " ) as h5_table:\n", - "\n", - " for data in create_stream(10):\n", - "\n", - " h5_table.write(\"table\", data)\n", - "\n", - " print(h5_table.h5file)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "So we can append some data groups. As long as the data group_name does not already exists. Let's try to overwrite the data group : data_1" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "try:\n", - " with HDF5TableWriter(\"container.h5\", mode=\"a\", group_name=\"data_1\") as h5_table:\n", - " for data in create_stream(10):\n", - " h5_table.write(\"table\", data)\n", - "except Exception as err:\n", - " print(\"Failed as expected:\", err)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Good ! I cannot overwrite my data." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(bool(h5_table.h5file.isopen))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Reading the Data" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Reading the whole table at once:\n", - "\n", - "For this, you have several choices. Since we used the HDF5TableWriter in this example, we have at least these options avilable:\n", - "\n", - "* Pandas\n", - "* PyTables\n", - "* Astropy Table\n", - "\n", - "For other TableWriter implementations, others may be possible (depending on format)\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Reading using `ctapipe.io.read_table`\n", - "\n", - "This is the preferred method, it returns an astropy `Table` and supports keeping track of\n", - "units, metadata and transformations." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.io import read_table\n", - "\n", - "\n", - "table = read_table(\"container.h5\", \"/data_0/table\")\n", - "table[:5]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "table.meta" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Reading with Pandas:\n", - "\n", - "Pandas is a convenient way to read the output. **HOWEVER BE WARNED** that so far Pandas does not support reading the table *meta-data* or *units* for colums, so that information is lost! " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import pandas as pd\n", - "\n", - "data = pd.read_hdf(\"container.h5\", key=\"/data_0/table\")\n", - "data.head()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Reading with PyTables" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import tables\n", - "\n", - "h5 = tables.open_file(\"container.h5\")\n", - "table = h5.root[\"data_0\"][\"table\"]\n", - "table" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "note that here we can still access the metadata" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "table.attrs" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Reading one-row-at-a-time:" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Rather than using the full-table methods, if you want to read it row-by-row (e.g. to maintain compatibility with an existing event loop), you can use a `TableReader` instance.\n", - "\n", - "The advantage here is that units and other metadata are retained and re-applied" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.io import HDF5TableReader\n", - "\n", - "\n", - "def read(mode):\n", - "\n", - " print(\"reading mode {}\".format(mode))\n", - "\n", - " with HDF5TableReader(\"container.h5\", mode=mode) as h5_table:\n", - "\n", - " for group_name in [\"data_0/\", \"data_1/\"]:\n", - "\n", - " group_name = \"/{}table\".format(group_name)\n", - " print(group_name)\n", - "\n", - " for data in h5_table.read(group_name, VariousTypesContainer):\n", - "\n", - " print(data.as_dict())" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "read(\"r\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "read(\"r+\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "read(\"a\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "read(\"w\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "language_info": { - "codemirror_mode": { - "name": "ipython" - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/docs/user-guide/examples_tutorials.rst b/docs/user-guide/examples_tutorials.rst new file mode 100644 index 00000000000..77b1161e010 --- /dev/null +++ b/docs/user-guide/examples_tutorials.rst @@ -0,0 +1,5 @@ +.. Include rendered examples build with sphinx-gallery +.. Actual files are in the examples/ directory in the +.. base of the repository + +.. include:: ../auto_examples/index.rst diff --git a/docs/user-guide/index.rst b/docs/user-guide/index.rst index 853408d8363..c74bf247684 100644 --- a/docs/user-guide/index.rst +++ b/docs/user-guide/index.rst @@ -9,6 +9,5 @@ User Guide getting-started tools data_models/index - tutorials/index - examples/index + examples_tutorials FAQ diff --git a/docs/user-guide/tutorials/calibrated_data_exploration.ipynb b/docs/user-guide/tutorials/calibrated_data_exploration.ipynb deleted file mode 100644 index 2966094d860..00000000000 --- a/docs/user-guide/tutorials/calibrated_data_exploration.ipynb +++ /dev/null @@ -1,398 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Explore Calibrated Data" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import ctapipe\n", - "from ctapipe.utils.datasets import get_dataset_path\n", - "from ctapipe.io import EventSource, EventSeeker\n", - "from ctapipe.visualization import CameraDisplay\n", - "from ctapipe.instrument import CameraGeometry\n", - "from matplotlib import pyplot as plt\n", - "from astropy import units as u\n", - "import numpy as np\n", - "\n", - "%matplotlib inline\n", - "plt.style.use(\"ggplot\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(ctapipe.__version__)\n", - "print(ctapipe.__file__)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's first open a raw event file and get an event out of it:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "filename = get_dataset_path(\"gamma_prod5.simtel.zst\")\n", - "source = EventSource(filename, max_events=2)\n", - "\n", - "for event in source:\n", - " print(event.index.event_id)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "filename" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "source" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "event" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "print(event.r1)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Perform basic calibration:\n", - "\n", - "Here we will use a `CameraCalibrator` which is just a simple wrapper that runs the three calibraraton and trace-integration phases of the pipeline, taking the data from levels:\n", - "\n", - " **R0** → **R1** → **DL0** → **DL1**\n", - "\n", - "You could of course do these each separately, by using the classes `R1Calibrator`, `DL0Reducer`, and `DL1Calibrator`.\n", - "Note that we have not specified any configuration to the `CameraCalibrator`, so it will be using the default algorithms and thresholds, other than specifying that the product is a \"HESSIOR1Calibrator\" (hopefully in the near future that will be automatic)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.calib import CameraCalibrator\n", - "\n", - "calib = CameraCalibrator(subarray=source.subarray)\n", - "calib(event)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now the *r1*, *dl0* and *dl1* containers are filled in the event\n", - "\n", - "* **r1.tel[x]**: contains the \"r1-calibrated\" waveforms, after gain-selection, pedestal subtraciton, and gain-correction\n", - "* **dl0.tel[x]**: is the same but with optional data volume reduction (some pixels not filled), in this case this is not performed by default, so it is the same as r1\n", - "* **dl1.tel[x]**: contains the (possibly re-calibrated) waveforms as dl0, but also the time-integrated *image* that has been calculated using a `ImageExtractor` (a `NeighborPeakWindowSum` by default)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for tel_id in event.dl1.tel:\n", - " print(\"TEL{:03}: {}\".format(tel_id, source.subarray.tel[tel_id]))\n", - " print(\" - r0 wave shape : {}\".format(event.r0.tel[tel_id].waveform.shape))\n", - " print(\" - r1 wave shape : {}\".format(event.r1.tel[tel_id].waveform.shape))\n", - " print(\" - dl1 image shape : {}\".format(event.dl1.tel[tel_id].image.shape))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Some image processing:\n", - "\n", - "Let's look at the image" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.visualization import CameraDisplay\n", - "\n", - "tel_id = sorted(event.r1.tel.keys())[1]\n", - "sub = source.subarray\n", - "geometry = sub.tel[tel_id].camera.geometry\n", - "image = event.dl1.tel[tel_id].image" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "disp = CameraDisplay(geometry, image=image)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.image import tailcuts_clean, hillas_parameters" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mask = tailcuts_clean(\n", - " geometry,\n", - " image,\n", - " picture_thresh=10,\n", - " boundary_thresh=5,\n", - " min_number_picture_neighbors=2,\n", - ")\n", - "cleaned = image.copy()\n", - "cleaned[~mask] = 0\n", - "disp = CameraDisplay(geometry, image=cleaned)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "params = hillas_parameters(geometry, cleaned)\n", - "print(params)\n", - "params" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "params = hillas_parameters(geometry, cleaned)\n", - "\n", - "plt.figure(figsize=(10, 10))\n", - "disp = CameraDisplay(geometry, image=image)\n", - "disp.add_colorbar()\n", - "disp.overlay_moments(params, color=\"red\", lw=3)\n", - "disp.highlight_pixels(mask, color=\"white\", alpha=0.3, linewidth=2)\n", - "\n", - "plt.xlim(params.x.to_value(u.m) - 0.5, params.x.to_value(u.m) + 0.5)\n", - "plt.ylim(params.y.to_value(u.m) - 0.5, params.y.to_value(u.m) + 0.5)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "source.metadata" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## More complex image processing:\n", - "\n", - "Let's now explore how stereo reconstruction works. \n", - "\n", - "### first, look at a summed image from multiple telescopes\n", - "\n", - "For this, we want to use a `CameraDisplay` again, but since we can't sum and display images with different cameras, we'll just sub-select images from a particular camera type\n", - "\n", - "These are the telescopes that are in this event:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tels_in_event = set(\n", - " event.dl1.tel.keys()\n", - ") # use a set here, so we can intersect it later\n", - "tels_in_event" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cam_ids = set(sub.get_tel_ids_for_type(\"MST_MST_NectarCam\"))\n", - "cam_ids" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cams_in_event = tels_in_event.intersection(cam_ids)\n", - "first_tel_id = list(cams_in_event)[0]\n", - "tel = sub.tel[first_tel_id]\n", - "print(\"{}s in event: {}\".format(tel, cams_in_event))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now, let's sum those images:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "image_sum = np.zeros_like(\n", - " tel.camera.geometry.pix_x.value\n", - ") # just make an array of 0's in the same shape as the camera\n", - "\n", - "for tel_id in cams_in_event:\n", - " image_sum += event.dl1.tel[tel_id].image" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "And finally display the sum of those images" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.figure(figsize=(8, 8))\n", - "\n", - "disp = CameraDisplay(tel.camera.geometry, image=image_sum)\n", - "disp.overlay_moments(params, with_label=False)\n", - "plt.title(\"Sum of {}x {}\".format(len(cams_in_event), tel))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "let's also show which telescopes those were. Note that currently ArrayDisplay's value field is a vector by `tel_index`, not `tel_id`, so we have to convert to a tel_index. (this may change in a future version to be more user-friendly)\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.visualization import ArrayDisplay" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "nectarcam_subarray = sub.select_subarray(cam_ids, name=\"NectarCam\")\n", - "\n", - "hit_pattern = np.zeros(shape=nectarcam_subarray.n_tels)\n", - "hit_pattern[[nectarcam_subarray.tel_indices[x] for x in cams_in_event]] = 100\n", - "\n", - "plt.set_cmap(plt.cm.Accent)\n", - "plt.figure(figsize=(8, 8))\n", - "\n", - "ad = ArrayDisplay(nectarcam_subarray)\n", - "ad.values = hit_pattern\n", - "ad.add_labels()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.8.13" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/docs/user-guide/tutorials/coordinates_example.ipynb b/docs/user-guide/tutorials/coordinates_example.ipynb deleted file mode 100644 index 5fa46d0850d..00000000000 --- a/docs/user-guide/tutorials/coordinates_example.ipynb +++ /dev/null @@ -1,649 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "# Coordinates usage in ctapipe" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "import astropy.units as u\n", - "import copy\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "\n", - "\n", - "from ctapipe.io import EventSource\n", - "from ctapipe.calib import CameraCalibrator\n", - "from ctapipe.utils import get_dataset_path\n", - "\n", - "from ctapipe.visualization import ArrayDisplay\n", - "\n", - "\n", - "%matplotlib inline" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "from astropy.coordinates import SkyCoord, AltAz\n", - "\n", - "from ctapipe.coordinates import (\n", - " GroundFrame,\n", - " TiltedGroundFrame,\n", - " NominalFrame,\n", - " TelescopeFrame,\n", - " CameraFrame,\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# make plots and fonts larger\n", - "plt.rcParams[\"figure.figsize\"] = (12, 8)\n", - "plt.rcParams[\"font.size\"] = 16" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "subslide" - } - }, - "source": [ - "## Open test dataset" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "filename = get_dataset_path(\"gamma_prod5.simtel.zst\")\n", - "source = EventSource(filename)\n", - "\n", - "events = [copy.deepcopy(event) for event in source]\n", - "event = events[4]\n", - "\n", - "layout = set(source.subarray.tel_ids)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Choose event with LST" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This ensures that the telescope is not \"parked\" (as it would be in an event where it is not triggered) but is actually pointing to a source." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "subslide" - } - }, - "outputs": [], - "source": [ - "print(f\"Telescope with data: {event.r1.tel.keys()}\")\n", - "tel_id = 3" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "## AltAz\n", - "\n", - "See [Astropy Docs on AltAz](https://docs.astropy.org/en/stable/api/astropy.coordinates.AltAz.html). \n", - "\n", - "Pointing direction of telescopes or the origin of a simulated shower are described in the `AltAz` frame.\n", - "This is a local, angular coordinate frame, with angles `altitude` and `azimuth`.\n", - "Altitude is the measured from the Horizon (0°) to the Zenith (90°).\n", - "For the azimuth, there are different conventions. In Astropy und thus ctapipe, Azimuth is oriented East of North (i.e., N=0°, E=90°)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from astropy.time import Time\n", - "from astropy.coordinates import EarthLocation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "obstime = Time(\"2013-11-01T03:00\")\n", - "location = EarthLocation.of_site(\"Roque de los Muchachos\")\n", - "\n", - "altaz = AltAz(location=location, obstime=obstime)\n", - "\n", - "array_pointing = SkyCoord(\n", - " alt=event.pointing.array_azimuth,\n", - " az=event.pointing.array_altitude,\n", - " frame=altaz,\n", - ")\n", - "\n", - "print(array_pointing)" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "subslide" - } - }, - "source": [ - "## CameraFrame\n", - "\n", - "Camera coordinate frame.\n", - "\n", - "The camera frame is a 2d cartesian frame, describing position of objects in the focal plane of the telescope.\n", - "\n", - "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.\n", - "\n", - "Now, x points north and y points west, so in this orientation, the camera coordinates line up with the CORSIKA ground coordinate system.\n", - "\n", - "MAGIC and FACT use a different camera coordinate system: Standing at the dish, looking at the camera, x points right, y points up.\n", - "To transform MAGIC/FACT to ctapipe, do x' = -y, y' = -x.\n", - "\n", - "**Typical usage**: Position of pixels in the focal plane." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "geometry = source.subarray.tel[tel_id].camera.geometry\n", - "pix_x = geometry.pix_x\n", - "pix_y = geometry.pix_y\n", - "focal_length = source.subarray.tel[tel_id].optics.equivalent_focal_length" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "telescope_pointing = SkyCoord(\n", - " alt=event.pointing.tel[tel_id].altitude,\n", - " az=event.pointing.tel[tel_id].azimuth,\n", - " frame=altaz,\n", - ")\n", - "\n", - "camera_frame = CameraFrame(\n", - " focal_length=focal_length,\n", - " rotation=0 * u.deg,\n", - " telescope_pointing=telescope_pointing,\n", - ")\n", - "\n", - "cam_coords = SkyCoord(x=pix_x, y=pix_y, frame=camera_frame)\n", - "\n", - "print(cam_coords)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "plt.scatter(cam_coords.x, cam_coords.y)\n", - "plt.title(f\"Camera type: {geometry.name}\")\n", - "plt.xlabel(f\"x / {cam_coords.x.unit}\")\n", - "plt.ylabel(f\"y / {cam_coords.y.unit}\")\n", - "plt.axis(\"square\");" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The implementation of the coordinate system with astropy makes it easier to use time of the observation and location of the observing site, to understand, for example which stars are visible during a certain night and how they might be visible in the camera.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.visualization import CameraDisplay\n", - "from ctapipe.instrument import SubarrayDescription\n", - "\n", - "location = EarthLocation.of_site(\"Roque de los Muchachos\")\n", - "obstime = Time(\"2018-11-01T04:00\")\n", - "\n", - "crab = SkyCoord.from_name(\"crab nebula\")\n", - "\n", - "altaz = AltAz(location=location, obstime=obstime)\n", - "\n", - "pointing = crab.transform_to(altaz)\n", - "\n", - "camera_frame = CameraFrame(\n", - " telescope_pointing=pointing,\n", - " focal_length=focal_length,\n", - " obstime=obstime,\n", - " location=location,\n", - ")\n", - "\n", - "\n", - "subarray = SubarrayDescription.read(\"dataset://gamma_prod5.simtel.zst\")\n", - "cam = subarray.tel[1].camera.geometry\n", - "fig, ax = plt.subplots()\n", - "display = CameraDisplay(cam, ax=ax)\n", - "\n", - "ax.set_title(\n", - " f\"La Palma, {obstime}, az={pointing.az.deg:.1f}°, zenith={pointing.zen.deg:.1f}°, camera={geometry.name}\"\n", - ")\n", - "\n", - "for i, name in enumerate([\"crab nebula\", \"o tau\", \"zet tau\"]):\n", - " star = SkyCoord.from_name(name)\n", - " star_cam = star.transform_to(camera_frame)\n", - "\n", - " x = star_cam.x.to_value(u.m)\n", - " y = star_cam.y.to_value(u.m)\n", - "\n", - " ax.plot(x, y, marker=\"*\", color=f\"C{i}\")\n", - " ax.annotate(\n", - " name,\n", - " xy=(x, y),\n", - " xytext=(5, 5),\n", - " textcoords=\"offset points\",\n", - " color=f\"C{i}\",\n", - " )\n", - "\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "subslide" - } - }, - "source": [ - "## TelescopeFrame\n", - "\n", - "Telescope coordinate frame.\n", - "A `Frame` using a `UnitSphericalRepresentation`.\n", - "\n", - "This is basically the same as a `HorizonCoordinate`, but the origin is at the telescope's pointing direction.\n", - "This is what astropy calls a `SkyOffsetFrame`.\n", - "\n", - "The axis of the telescope frame, `fov_lon` and `fov_lat`, are aligned with the horizontal system's azimuth and altitude respectively.\n", - " \n", - "Pointing corrections should applied to the transformation between this frame and the camera frame." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true, - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "telescope_frame = TelescopeFrame(\n", - " telescope_pointing=pointing,\n", - " obstime=pointing.obstime,\n", - " location=pointing.location,\n", - ")\n", - "telescope_coords = cam_coords.transform_to(telescope_frame)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "wrap_angle = telescope_pointing.az + 180 * u.deg\n", - "\n", - "plt.axis(\"equal\")\n", - "plt.scatter(\n", - " telescope_coords.fov_lon.deg, telescope_coords.fov_lat.deg, alpha=0.2, color=\"gray\"\n", - ")\n", - "\n", - "\n", - "for i, name in enumerate([\"crab nebula\", \"o tau\", \"zet tau\"]):\n", - " star = SkyCoord.from_name(name)\n", - " star_tel = star.transform_to(telescope_frame)\n", - "\n", - " plt.plot(star_tel.fov_lon.deg, star_tel.fov_lat.deg, \"*\", ms=10)\n", - " plt.annotate(\n", - " name,\n", - " xy=(star_tel.fov_lon.deg, star_tel.fov_lat.deg),\n", - " xytext=(5, 5),\n", - " textcoords=\"offset points\",\n", - " color=f\"C{i}\",\n", - " )\n", - "\n", - "plt.xlabel(\"fov_lon / {}\".format(telescope_coords.altaz.az.unit))\n", - "plt.ylabel(\"fov_lat / {}\".format(telescope_coords.altaz.alt.unit))" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "subslide" - } - }, - "source": [ - "## NominalFrame\n", - "\n", - "Nominal coordinate frame.\n", - "A Frame using a `UnitSphericalRepresentation`.\n", - "This is basically the same as a `HorizonCoordinate`, but the\n", - "origin is at an arbitray position in the sky.\n", - "This is what astropy calls a `SkyOffsetFrame`\n", - "If the telescopes are in divergent pointing, this `Frame` can be\n", - "used to transform to a common system.\n", - "- 2D reconstruction (`HillasIntersector`) is performed in this frame \n", - "- 3D reconstruction (`HillasReconstructor`) doesn't need this frame" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "source": [ - "Let's play a bit with 3 LSTs with divergent pointing" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "location = EarthLocation.of_site(\"Roque de los Muchachos\")\n", - "obstime = Time(\"2018-11-01T02:00\")\n", - "altaz = AltAz(location=location, obstime=obstime)\n", - "\n", - "crab = SkyCoord.from_name(\"crab nebula\")\n", - "\n", - "# let's observe crab\n", - "array_pointing = crab.transform_to(altaz)\n", - "\n", - "\n", - "# let the telescopes point to different positions\n", - "alt_offsets = u.Quantity([1, -1, -1], u.deg)\n", - "az_offsets = u.Quantity([0, -2, +2], u.deg)\n", - "\n", - "\n", - "tel_pointings = SkyCoord(\n", - " alt=array_pointing.alt + alt_offsets,\n", - " az=array_pointing.az + az_offsets,\n", - " frame=altaz,\n", - ")\n", - "\n", - "camera_frames = CameraFrame(\n", - " telescope_pointing=tel_pointings, # multiple pointings, so we get multiple frames\n", - " focal_length=focal_length,\n", - " obstime=obstime,\n", - " location=location,\n", - ")\n", - "\n", - "nom_frame = NominalFrame(origin=array_pointing, obstime=obstime, location=location)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(15, 10))\n", - "ax.set_aspect(1)\n", - "\n", - "for i in range(3):\n", - " cam_coord = SkyCoord(x=pix_x, y=pix_y, frame=camera_frames[i])\n", - " nom_coord = cam_coord.transform_to(nom_frame)\n", - "\n", - " ax.scatter(\n", - " x=nom_coord.fov_lon.deg,\n", - " y=nom_coord.fov_lat.deg,\n", - " label=f\"Telescope {i + 1}\",\n", - " s=30,\n", - " alpha=0.15,\n", - " )\n", - "\n", - "\n", - "for i, name in enumerate([\"Crab\", \"o tau\", \"zet tau\"]):\n", - " s = SkyCoord.from_name(name)\n", - " s_nom = s.transform_to(nom_frame)\n", - " ax.plot(\n", - " s_nom.fov_lon.deg,\n", - " s_nom.fov_lat.deg,\n", - " \"*\",\n", - " ms=10,\n", - " )\n", - " ax.annotate(\n", - " name,\n", - " xy=(s_nom.fov_lon.deg, s_nom.fov_lat.deg),\n", - " xytext=(5, 5),\n", - " textcoords=\"offset points\",\n", - " color=f\"C{i}\",\n", - " )\n", - "\n", - "\n", - "ax.set_xlabel(f\"fov_lon / deg\")\n", - "ax.set_ylabel(f\"fov_lat / deg\")\n", - "\n", - "ax.legend()\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "subslide" - } - }, - "source": [ - "## GroundFrame\n", - "\n", - "\n", - "Ground coordinate frame. The ground coordinate frame is a simple\n", - " cartesian frame describing the 3 dimensional position of objects\n", - " compared to the array ground level in relation to the nomial\n", - " centre of the array. Typically this frame will be used for\n", - " describing the position on telescopes and equipment\n", - " \n", - "**Typical usage**: positions of telescopes on the ground (x, y, z)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "source.subarray.peek()" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "source": [ - "In case a layout is selected, the following line will produce a different output from the picture above." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "source.subarray.select_subarray(layout, name=\"Prod3b layout\").peek()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "![Ground Frame](ground_frame.png)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In this image all the telescope from the `gamma_test.simtel.gz` file are plotted as spheres in the GroundFrame." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## TiltedGroundFrame" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Tilted ground coordinate frame. \n", - "\n", - "The tilted ground coordinate frame is a cartesian system describing the 2 dimensional projected positions of objects in a tilted plane described by pointing_direction. The plane is rotated along the z_axis by the azimuth of the `pointing_direction` and then it is inclined with an angle equal to the zenith angle of the `pointing_direction`.\n", - "\n", - "This frame is used for the reconstruction of the shower core position." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "![Tilted Ground Frame](tilted_ground_frame.png)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This image picture both the telescopes in the GroundFrame (red) and in the TiltedGroundFrame (green) are displayed: in this case since the azimuth of the `pointing_direction` is 0 degrees, then the plane is just tilted according to the zenith angle. \n", - "\n", - "For playing with these and with more 3D models of the telescopes themselves, have a look at the [CREED_VTK](https://github.com/thomasgas/CREED_VTK) library. " - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.13" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/docs/user-guide/tutorials/ctapipe_handson.ipynb b/docs/user-guide/tutorials/ctapipe_handson.ipynb deleted file mode 100644 index f3c27873309..00000000000 --- a/docs/user-guide/tutorials/ctapipe_handson.ipynb +++ /dev/null @@ -1,678 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Getting Started with ctapipe\n", - "\n", - "This hands-on was presented at the Paris CTA Consoritum meeting (K. Kosack)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Part 1: load and loop over data" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.io import EventSource\n", - "from ctapipe import utils\n", - "from matplotlib import pyplot as plt\n", - "import numpy as np\n", - "\n", - "%matplotlib inline" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "path = utils.get_dataset_path(\"gamma_prod5.simtel.zst\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "source = EventSource(path, max_events=5)\n", - "\n", - "for event in source:\n", - " print(event.count, event.index.event_id, event.simulation.shower.energy)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "event" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "event.r1" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for event in EventSource(path, max_events=5):\n", - " print(event.count, event.r1.tel.keys())" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "event.r0.tel[3]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "r0tel = event.r0.tel[3]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "r0tel.waveform" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "r0tel.waveform.shape" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "note that this is ($N_{channels}$, $N_{pixels}$, $N_{samples}$)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.pcolormesh(r0tel.waveform[0])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "brightest_pixel = np.argmax(r0tel.waveform[0].sum(axis=1))\n", - "print(f\"pixel {brightest_pixel} has sum {r0tel.waveform[0,1535].sum()}\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.plot(r0tel.waveform[0, brightest_pixel], label=\"channel 0 (high-gain)\")\n", - "plt.plot(r0tel.waveform[1, brightest_pixel], label=\"channel 1 (low-gain)\")\n", - "plt.legend()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ipywidgets import interact\n", - "\n", - "\n", - "@interact\n", - "def view_waveform(chan=0, pix_id=brightest_pixel):\n", - " plt.plot(r0tel.waveform[chan, pix_id])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "try making this compare 2 waveforms" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Part 2: Explore the instrument description\n", - "This is all well and good, but we don't really know what camera or telescope this is... how do we get instrumental description info?\n", - "\n", - "Currently this is returned *inside* the event (it will soon change to be separate in next version or so)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray = source.subarray" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.peek()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.to_table()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.tel[2]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.tel[2].camera" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.tel[2].optics" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tel = subarray.tel[2]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tel.camera" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "tel.optics" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tel.camera.geometry.pix_x" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tel.camera.geometry.to_table()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tel.optics.mirror_area" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.visualization import CameraDisplay" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "disp = CameraDisplay(tel.camera.geometry)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "disp = CameraDisplay(tel.camera.geometry)\n", - "disp.image = r0tel.waveform[\n", - " 0, :, 10\n", - "] # display channel 0, sample 0 (try others like 10)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - " ** aside: ** show demo using a CameraDisplay in interactive mode in ipython rather than notebook" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Part 3: Apply some calibration and trace integration" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.calib import CameraCalibrator" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "calib = CameraCalibrator(subarray=subarray)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for event in EventSource(path, max_events=5):\n", - " calib(event) # fills in r1, dl0, and dl1\n", - " print(event.dl1.tel.keys())" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "event.dl1.tel[3]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "dl1tel = event.dl1.tel[3]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "dl1tel.image.shape # note this will be gain-selected in next version, so will be just 1D array of 1855" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "dl1tel.peak_time" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "CameraDisplay(tel.camera.geometry, image=dl1tel.image)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "CameraDisplay(tel.camera.geometry, image=dl1tel.peak_time)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now for Hillas Parameters" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.image import hillas_parameters, tailcuts_clean" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "image = dl1tel.image\n", - "mask = tailcuts_clean(tel.camera.geometry, image, picture_thresh=10, boundary_thresh=5)\n", - "mask" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "CameraDisplay(tel.camera.geometry, image=mask)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cleaned = image.copy()\n", - "cleaned[~mask] = 0" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "disp = CameraDisplay(tel.camera.geometry, image=cleaned)\n", - "disp.cmap = plt.cm.coolwarm\n", - "disp.add_colorbar()\n", - "plt.xlim(0.5, 1.0)\n", - "plt.ylim(-1.0, 0.0)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "params = hillas_parameters(tel.camera.geometry, cleaned)\n", - "print(params)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "disp = CameraDisplay(tel.camera.geometry, image=cleaned)\n", - "disp.cmap = plt.cm.coolwarm\n", - "disp.add_colorbar()\n", - "plt.xlim(0.5, 1.0)\n", - "plt.ylim(-1.0, 0.0)\n", - "disp.overlay_moments(params, color=\"white\", lw=2)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Part 4: Let's put it all together: \n", - "- loop over events, selecting only telescopes of the same type (e.g. LST:LSTCam)\n", - "- for each event, apply calibration/trace integration\n", - "- calculate Hillas parameters \n", - "- write out all hillas paremeters to a file that can be loaded with Pandas" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "first let's select only those telescopes with LST:LSTCam" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.telescope_types" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "subarray.get_tel_ids_for_type(\"LST_LST_LSTCam\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now let's write out program" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "data = utils.get_dataset_path(\"gamma_prod5.simtel.zst\")\n", - "source = EventSource(data) # remove the max_events limit to get more stats" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for event in source:\n", - " calib(event)\n", - "\n", - " for tel_id, tel_data in event.dl1.tel.items():\n", - " tel = source.subarray.tel[tel_id]\n", - " mask = tailcuts_clean(tel.camera.geometry, tel_data.image)\n", - " if np.count_nonzero(mask) > 0:\n", - " params = hillas_parameters(tel.camera.geometry[mask], tel_data.image[mask])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.io import HDF5TableWriter" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "with HDF5TableWriter(filename=\"hillas.h5\", group_name=\"dl1\", overwrite=True) as writer:\n", - "\n", - " source = EventSource(data, allowed_tels=[1, 2, 3, 4], max_events=10)\n", - " for event in source:\n", - " calib(event)\n", - "\n", - " for tel_id, tel_data in event.dl1.tel.items():\n", - " tel = source.subarray.tel[tel_id]\n", - " mask = tailcuts_clean(tel.camera.geometry, tel_data.image)\n", - " params = hillas_parameters(tel.camera.geometry[mask], tel_data.image[mask])\n", - " writer.write(\"hillas\", params)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### We can now load in the file we created and plot it" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "!ls *.h5" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import pandas as pd\n", - "\n", - "hillas = pd.read_hdf(\"hillas.h5\", key=\"/dl1/hillas\")\n", - "hillas" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "_ = hillas.hist(figsize=(8, 8))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "If you do this yourself, chose a larger file to loop over more events to get better statistics" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.8.13" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/docs/user-guide/tutorials/ctapipe_overview.ipynb b/docs/user-guide/tutorials/ctapipe_overview.ipynb deleted file mode 100644 index 6ee1bac2b63..00000000000 --- a/docs/user-guide/tutorials/ctapipe_overview.ipynb +++ /dev/null @@ -1,1180 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": { - "pycharm": {} - }, - "source": [ - "# Analyzing Events Using ctapipe" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {}, - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "
\n", - "\n", - "\"ctapipe\"/\n", - "\n", - "\n", - "

Initially presented @ LST Analysis Bootcamp

\n", - "\n", - "

Padova, 26.11.2018

\n", - "\n", - "

Maximilian Nöthe (@maxnoe) & Kai A. Brügge (@mackaiver)

\n", - "\n", - "
" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {}, - "slideshow": { - "slide_type": "skip" - } - }, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "\n", - "%matplotlib inline" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {}, - "slideshow": { - "slide_type": "skip" - } - }, - "outputs": [], - "source": [ - "plt.rcParams[\"figure.figsize\"] = (12, 8)\n", - "plt.rcParams[\"font.size\"] = 14\n", - "plt.rcParams[\"figure.figsize\"]" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {}, - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "

Table of Contents

\n", - "
" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {}, - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "## General Information" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {}, - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "### Design\n", - "\n", - "* DL0 → DL3 analysis\n", - "\n", - "* Currently some R0 → DL2 code to be able to analyze simtel files\n", - "\n", - "* ctapipe is built upon the Scientific Python Stack, core dependencies are\n", - " * numpy\n", - " * scipy\n", - " * astropy\n", - " * numba" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {}, - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "### Developement\n", - "\n", - "* ctapipe is developed as Open Source Software (BSD 3-Clause License) at \n", - "\n", - "* We use the \"Github-Workflow\": \n", - " * Few people (e.g. @kosack, @maxnoe) have write access to the main repository\n", - " * Contributors fork the main repository and work on branches\n", - " * Pull Requests are merged after Code Review and automatic execution of the test suite\n", - "\n", - "* Early developement stage ⇒ backwards-incompatible API changes might and will happen \n" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {}, - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "### What's there?\n", - "\n", - "* Reading simtel simulation files\n", - "* Simple calibration, cleaning and feature extraction functions\n", - "* Camera and Array plotting\n", - "* Coordinate frames and transformations \n", - "* Stereo-reconstruction using line intersections\n", - " \n", - " " - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {}, - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "### What's still missing?\n", - "\n", - "* Good integration with machine learning techniques\n", - "* IRF calculation \n", - "* Documentation, e.g. formal definitions of coordinate frames " - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {}, - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "### What can you do?\n", - "\n", - "* Report issues\n", - " * Hard to get started? Tell us where you are stuck\n", - " * Tell user stories\n", - " * Missing features\n", - "\n", - "* Start contributing\n", - " * ctapipe needs more workpower\n", - " * Implement new reconstruction features" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {}, - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "## A simple hillas analysis" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {} - }, - "source": [ - "### Reading in simtel files" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "from ctapipe.io import EventSource\n", - "from ctapipe.utils.datasets import get_dataset_path\n", - "\n", - "input_url = get_dataset_path(\"gamma_prod5.simtel.zst\")\n", - "\n", - "# EventSource() automatically detects what kind of file we are giving it,\n", - "# if already supported by ctapipe\n", - "source = EventSource(input_url, max_events=5)\n", - "\n", - "print(type(source))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "for event in source:\n", - " print(\n", - " \"Id: {}, E = {:1.3f}, Telescopes: {}\".format(\n", - " event.count, event.simulation.shower.energy, len(event.r0.tel)\n", - " )\n", - " )" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {} - }, - "source": [ - "Each event is a `DataContainer` holding several `Field`s of data, which can be containers or just numbers.\n", - "Let's look a one event:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "event" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "source.subarray.camera_types" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "len(event.r0.tel), len(event.r1.tel)" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {} - }, - "source": [ - "### Data calibration\n", - "\n", - "The `CameraCalibrator` calibrates the event (obtaining the `dl1` images)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "from ctapipe.calib import CameraCalibrator\n", - "\n", - "calibrator = CameraCalibrator(subarray=source.subarray)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "calibrator(event)" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {} - }, - "source": [ - "### Event displays\n", - "\n", - "Let's use ctapipe's plotting facilities to plot the telescope images" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "event.dl1.tel.keys()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "tel_id = 130" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "geometry = source.subarray.tel[tel_id].camera.geometry\n", - "dl1 = event.dl1.tel[tel_id]\n", - "\n", - "geometry, dl1" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "dl1.image" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "from ctapipe.visualization import CameraDisplay\n", - "\n", - "display = CameraDisplay(geometry)\n", - "\n", - "# right now, there might be one image per gain channel.\n", - "# This will change as soon as\n", - "display.image = dl1.image\n", - "display.add_colorbar()" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {} - }, - "source": [ - "### Image Cleaning" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "from ctapipe.image.cleaning import tailcuts_clean" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "# unoptimized cleaning levels\n", - "cleaning_level = {\n", - " \"CHEC\": (2, 4, 2),\n", - " \"LSTCam\": (3.5, 7, 2),\n", - " \"FlashCam\": (3.5, 7, 2),\n", - " \"NectarCam\": (4, 8, 2),\n", - "}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "boundary, picture, min_neighbors = cleaning_level[geometry.name]\n", - "\n", - "clean = tailcuts_clean(\n", - " geometry,\n", - " dl1.image,\n", - " boundary_thresh=boundary,\n", - " picture_thresh=picture,\n", - " min_number_picture_neighbors=min_neighbors,\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))\n", - "\n", - "d1 = CameraDisplay(geometry, ax=ax1)\n", - "d2 = CameraDisplay(geometry, ax=ax2)\n", - "\n", - "ax1.set_title(\"Image\")\n", - "d1.image = dl1.image\n", - "d1.add_colorbar(ax=ax1)\n", - "\n", - "ax2.set_title(\"Pulse Time\")\n", - "d2.image = dl1.peak_time - np.average(dl1.peak_time, weights=dl1.image)\n", - "d2.cmap = \"RdBu_r\"\n", - "d2.add_colorbar(ax=ax2)\n", - "d2.set_limits_minmax(-20, 20)\n", - "\n", - "d1.highlight_pixels(clean, color=\"red\", linewidth=1)" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {} - }, - "source": [ - "### Image Parameters" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "from ctapipe.image import (\n", - " hillas_parameters,\n", - " leakage_parameters,\n", - " concentration_parameters,\n", - ")\n", - "from ctapipe.image import timing_parameters\n", - "from ctapipe.image import number_of_islands\n", - "from ctapipe.image import camera_to_shower_coordinates" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "hillas = hillas_parameters(geometry[clean], dl1.image[clean])\n", - "\n", - "print(hillas)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "display = CameraDisplay(geometry)\n", - "\n", - "# set \"unclean\" pixels to 0\n", - "cleaned = dl1.image.copy()\n", - "cleaned[~clean] = 0.0\n", - "\n", - "display.image = cleaned\n", - "display.add_colorbar()\n", - "\n", - "display.overlay_moments(hillas, color=\"xkcd:red\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "timing = timing_parameters(geometry, dl1.image, dl1.peak_time, hillas, clean)\n", - "\n", - "print(timing)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "long, trans = camera_to_shower_coordinates(\n", - " geometry.pix_x, geometry.pix_y, hillas.x, hillas.y, hillas.psi\n", - ")\n", - "\n", - "plt.plot(long[clean], dl1.peak_time[clean], \"o\")\n", - "plt.plot(long[clean], timing.slope * long[clean] + timing.intercept)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "l = leakage_parameters(geometry, dl1.image, clean)\n", - "print(l)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "disp = CameraDisplay(geometry)\n", - "disp.image = dl1.image\n", - "disp.highlight_pixels(geometry.get_border_pixel_mask(1), linewidth=2, color=\"xkcd:red\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "n_islands, island_id = number_of_islands(geometry, clean)\n", - "\n", - "print(n_islands)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "conc = concentration_parameters(geometry, dl1.image, hillas)\n", - "print(conc)" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {} - }, - "source": [ - "### Putting it all together / Stereo reconstruction\n", - "\n", - "\n", - "All these steps are now unified in several components configurable through the config system, mainly:\n", - "\n", - "* CameraCalibrator for DL0 → DL1 (Images)\n", - "* ImageProcessor for DL1 (Images) → DL1 (Parameters)\n", - "* ShowerProcessor for stereo reconstruction of the shower geometry\n", - "* DataWriter for writing data into HDF5\n", - "\n", - "A command line tool doing these steps and writing out data in HDF5 format is available as `ctapipe-process`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {}, - "scrolled": false - }, - "outputs": [], - "source": [ - "import astropy.units as u\n", - "from astropy.coordinates import SkyCoord, AltAz\n", - "\n", - "from ctapipe.containers import ImageParametersContainer\n", - "from ctapipe.io import EventSource\n", - "from ctapipe.utils.datasets import get_dataset_path\n", - "\n", - "from ctapipe.calib import CameraCalibrator\n", - "\n", - "from ctapipe.image import ImageProcessor\n", - "from ctapipe.reco import ShowerProcessor\n", - "\n", - "from ctapipe.io import DataWriter\n", - "\n", - "from copy import deepcopy\n", - "import tempfile\n", - "\n", - "from traitlets.config import Config\n", - "\n", - "\n", - "image_processor_config = Config(\n", - " {\n", - " \"ImageProcessor\": {\n", - " \"image_cleaner_type\": \"TailcutsImageCleaner\",\n", - " \"TailcutsImageCleaner\": {\n", - " \"picture_threshold_pe\": [\n", - " (\"type\", \"LST_LST_LSTCam\", 7.5),\n", - " (\"type\", \"MST_MST_FlashCam\", 8),\n", - " (\"type\", \"MST_MST_NectarCam\", 8),\n", - " (\"type\", \"SST_ASTRI_CHEC\", 7),\n", - " ],\n", - " \"boundary_threshold_pe\": [\n", - " (\"type\", \"LST_LST_LSTCam\", 5),\n", - " (\"type\", \"MST_MST_FlashCam\", 4),\n", - " (\"type\", \"MST_MST_NectarCam\", 4),\n", - " (\"type\", \"SST_ASTRI_CHEC\", 4),\n", - " ],\n", - " },\n", - " }\n", - " }\n", - ")\n", - "\n", - "input_url = get_dataset_path(\"gamma_prod5.simtel.zst\")\n", - "source = EventSource(input_url)\n", - "\n", - "calibrator = CameraCalibrator(subarray=source.subarray)\n", - "image_processor = ImageProcessor(\n", - " subarray=source.subarray, config=image_processor_config\n", - ")\n", - "shower_processor = ShowerProcessor(subarray=source.subarray)\n", - "horizon_frame = AltAz()\n", - "\n", - "f = tempfile.NamedTemporaryFile(suffix=\".hdf5\")\n", - "\n", - "with DataWriter(\n", - " source, output_path=f.name, overwrite=True, write_showers=True\n", - ") as writer:\n", - "\n", - " for event in source:\n", - " energy = event.simulation.shower.energy\n", - " n_telescopes_r1 = len(event.r1.tel)\n", - " event_id = event.index.event_id\n", - " print(f\"Id: {event_id}, E = {energy:1.3f}, Telescopes (R1): {n_telescopes_r1}\")\n", - "\n", - " calibrator(event)\n", - " image_processor(event)\n", - " shower_processor(event)\n", - "\n", - " stereo = event.dl2.stereo.geometry[\"HillasReconstructor\"]\n", - " if stereo.is_valid:\n", - " print(\" Alt: {:.2f}°\".format(stereo.alt.deg))\n", - " print(\" Az: {:.2f}°\".format(stereo.az.deg))\n", - " print(\" Hmax: {:.0f}\".format(stereo.h_max))\n", - " print(\" CoreX: {:.1f}\".format(stereo.core_x))\n", - " print(\" CoreY: {:.1f}\".format(stereo.core_y))\n", - " print(\" Multiplicity: {:d}\".format(len(stereo.telescopes)))\n", - "\n", - " # save a nice event for plotting later\n", - " if event.count == 3:\n", - " plotting_event = deepcopy(event)\n", - "\n", - " writer(event)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "from astropy.coordinates.angle_utilities import angular_separation\n", - "import pandas as pd\n", - "\n", - "from ctapipe.io import TableLoader\n", - "\n", - "loader = TableLoader(f.name, load_dl2=True, load_simulated=True)\n", - "\n", - "events = loader.read_subarray_events()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "theta = angular_separation(\n", - " events[\"HillasReconstructor_az\"].quantity,\n", - " events[\"HillasReconstructor_alt\"].quantity,\n", - " events[\"true_az\"].quantity,\n", - " events[\"true_alt\"].quantity,\n", - ")\n", - "\n", - "plt.hist(theta.to_value(u.deg) ** 2, bins=25, range=[0, 0.3])\n", - "plt.xlabel(r\"$\\theta² / deg²$\")\n", - "None" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {} - }, - "source": [ - "## ArrayDisplay\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "from ctapipe.visualization import ArrayDisplay\n", - "\n", - "\n", - "angle_offset = plotting_event.pointing.array_azimuth\n", - "\n", - "plotting_hillas = {\n", - " tel_id: dl1.parameters.hillas for tel_id, dl1 in plotting_event.dl1.tel.items()\n", - "}\n", - "\n", - "plotting_core = {\n", - " tel_id: dl1.parameters.core.psi for tel_id, dl1 in plotting_event.dl1.tel.items()\n", - "}\n", - "\n", - "\n", - "disp = ArrayDisplay(source.subarray)\n", - "\n", - "disp.set_line_hillas(plotting_hillas, plotting_core, 500)\n", - "\n", - "plt.scatter(\n", - " plotting_event.simulation.shower.core_x,\n", - " plotting_event.simulation.shower.core_y,\n", - " s=200,\n", - " c=\"k\",\n", - " marker=\"x\",\n", - " label=\"True Impact\",\n", - ")\n", - "plt.scatter(\n", - " plotting_event.dl2.stereo.geometry[\"HillasReconstructor\"].core_x,\n", - " plotting_event.dl2.stereo.geometry[\"HillasReconstructor\"].core_y,\n", - " s=200,\n", - " c=\"r\",\n", - " marker=\"x\",\n", - " label=\"Estimated Impact\",\n", - ")\n", - "\n", - "plt.legend()\n", - "# plt.xlim(-400, 400)\n", - "# plt.ylim(-400, 400)" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {} - }, - "source": [ - "### Reading the LST dl1 data\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "loader = TableLoader(f.name, load_simulated=True, load_dl1_parameters=True)\n", - "\n", - "dl1_table = loader.read_telescope_events([\"LST_LST_LSTCam\"])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "plt.scatter(\n", - " np.log10(dl1_table[\"true_energy\"].quantity / u.TeV),\n", - " np.log10(dl1_table[\"hillas_intensity\"]),\n", - ")\n", - "plt.xlabel(\"log10(E / TeV)\")\n", - "plt.ylabel(\"log10(intensity)\")\n", - "None" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {} - }, - "source": [ - "## Isn't python slow?\n", - "\n", - "* Many of you might have heard: \"Python is slow\".\n", - "* That's trueish.\n", - "* All python objects are classes living on the heap, even integers.\n", - "* Looping over lots of \"primitives\" is quite slow compared to other languages.\n", - "\n", - "⇒ Vectorize as much as possible using numpy \n", - "⇒ Use existing interfaces to fast C / C++ / Fortran code \n", - "⇒ Optimize using numba \n", - "\n", - "**But: \"Premature Optimization is the root of all evil\" — Donald Knuth**\n", - "\n", - "So profile to find exactly what is slow.\n", - "\n", - "### Why use python then?\n", - "\n", - "* Python works very well as *glue* for libraries of all kinds of languages\n", - "* Python has a rich ecosystem for data science, physics, algorithms, astronomy\n", - "\n", - "### Example: Number of Islands\n", - "\n", - "Find all groups of pixels, that survived the cleaning" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "from ctapipe.image import toymodel\n", - "from ctapipe.instrument import SubarrayDescription\n", - "\n", - "\n", - "geometry = loader.subarray.tel[1].camera.geometry" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {} - }, - "source": [ - "Let's create a toy images with several islands;" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "np.random.seed(42)\n", - "\n", - "image = np.zeros(geometry.n_pixels)\n", - "\n", - "\n", - "for i in range(9):\n", - "\n", - " model = toymodel.Gaussian(\n", - " x=np.random.uniform(-0.8, 0.8) * u.m,\n", - " y=np.random.uniform(-0.8, 0.8) * u.m,\n", - " width=np.random.uniform(0.05, 0.075) * u.m,\n", - " length=np.random.uniform(0.1, 0.15) * u.m,\n", - " psi=np.random.uniform(0, 2 * np.pi) * u.rad,\n", - " )\n", - "\n", - " new_image, sig, bg = model.generate_image(\n", - " geometry, intensity=np.random.uniform(1000, 3000), nsb_level_pe=5\n", - " )\n", - " image += new_image" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "clean = tailcuts_clean(\n", - " geometry,\n", - " image,\n", - " picture_thresh=10,\n", - " boundary_thresh=5,\n", - " min_number_picture_neighbors=2,\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "disp = CameraDisplay(geometry)\n", - "disp.image = image\n", - "disp.highlight_pixels(clean, color=\"xkcd:red\", linewidth=1.5)\n", - "disp.add_colorbar()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "def num_islands_python(camera, clean):\n", - " \"\"\"A breadth first search to find connected islands of neighboring pixels in the cleaning set\"\"\"\n", - "\n", - " # the camera geometry has a [n_pixel, n_pixel] boolean array\n", - " # that is True where two pixels are neighbors\n", - " neighbors = camera.neighbor_matrix\n", - "\n", - " island_ids = np.zeros(camera.n_pixels)\n", - " current_island = 0\n", - "\n", - " # a set to remember which pixels we already visited\n", - " visited = set()\n", - "\n", - " # go only through the pixels, that survived cleaning\n", - " for pix_id in np.where(clean)[0]:\n", - " if pix_id not in visited:\n", - " # remember that we already checked this pixel\n", - " visited.add(pix_id)\n", - "\n", - " # if we land in the outer loop again, we found a new island\n", - " current_island += 1\n", - " island_ids[pix_id] = current_island\n", - "\n", - " # now check all neighbors of the current pixel recursively\n", - " to_check = set(np.where(neighbors[pix_id] & clean)[0])\n", - " while to_check:\n", - " pix_id = to_check.pop()\n", - "\n", - " if pix_id not in visited:\n", - " visited.add(pix_id)\n", - " island_ids[pix_id] = current_island\n", - "\n", - " to_check.update(np.where(neighbors[pix_id] & clean)[0])\n", - "\n", - " n_islands = current_island\n", - " return n_islands, island_ids" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "n_islands, island_ids = num_islands_python(geometry, clean)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "from matplotlib.colors import ListedColormap\n", - "\n", - "cmap = plt.get_cmap(\"Paired\")\n", - "cmap = ListedColormap(cmap.colors[:n_islands])\n", - "cmap.set_under(\"k\")\n", - "\n", - "disp = CameraDisplay(geometry)\n", - "disp.image = island_ids\n", - "disp.cmap = cmap\n", - "disp.set_limits_minmax(0.5, n_islands + 0.5)\n", - "disp.add_colorbar()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "%timeit num_islands_python(geometry, clean)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "from scipy.sparse.csgraph import connected_components\n", - "\n", - "\n", - "def num_islands_scipy(geometry, clean):\n", - " neighbors = geometry.neighbor_matrix_sparse\n", - "\n", - " clean_neighbors = neighbors[clean][:, clean]\n", - " num_islands, labels = connected_components(clean_neighbors, directed=False)\n", - "\n", - " island_ids = np.zeros(geometry.n_pixels)\n", - " island_ids[clean] = labels + 1\n", - "\n", - " return num_islands, island_ids" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "n_islands_s, island_ids_s = num_islands_scipy(geometry, clean)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "disp = CameraDisplay(geometry)\n", - "disp.image = island_ids_s\n", - "disp.cmap = cmap\n", - "disp.set_limits_minmax(0.5, n_islands_s + 0.5)\n", - "disp.add_colorbar()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "pycharm": {} - }, - "outputs": [], - "source": [ - "%timeit num_islands_scipy(geometry, clean)" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "pycharm": {} - }, - "source": [ - "**A lot less code, and a factor 3 speed improvement**" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, current ctapipe implementation is using numba:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%timeit number_of_islands(geometry, clean)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.13" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/docs/user-guide/tutorials/index.rst b/docs/user-guide/tutorials/index.rst deleted file mode 100644 index bfa4bb78594..00000000000 --- a/docs/user-guide/tutorials/index.rst +++ /dev/null @@ -1,14 +0,0 @@ -.. _tutorials: - -========= -Tutorials -========= - -.. toctree:: - - ctapipe_handson - ctapipe_overview - coordinates_example - raw_data_exploration - calibrated_data_exploration - theta_square diff --git a/docs/user-guide/tutorials/raw_data_exploration.ipynb b/docs/user-guide/tutorials/raw_data_exploration.ipynb deleted file mode 100644 index 1f2a23c978b..00000000000 --- a/docs/user-guide/tutorials/raw_data_exploration.ipynb +++ /dev/null @@ -1,522 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Exploring Raw Data\n", - "\n", - "Here are just some very simple examples of going through and inspecting the raw data, and making some plots using `ctapipe`.\n", - "The data explored here are *raw Monte Carlo* data, which is Data Level \"R0\" in CTA terminology (e.g. it is before any processing that would happen inside a Camera or off-line)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Setup:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from ctapipe.utils import get_dataset_path\n", - "from ctapipe.io import EventSource\n", - "from ctapipe.visualization import CameraDisplay\n", - "from ctapipe.instrument import CameraGeometry\n", - "from matplotlib import pyplot as plt\n", - "from astropy import units as u\n", - "\n", - "%matplotlib inline" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To read SimTelArray format data, ctapipe uses the `pyeventio` library (which is installed automatically along with ctapipe). The following lines however will load any data known to ctapipe (multiple `EventSources` are implemented, and chosen automatically based on the type of the input file. \n", - "\n", - "All data access first starts with an `EventSource`, and here we use a helper function `event_source` that constructs one. The resulting `source` object can be iterated over like a list of events. We also here use an `EventSeeker` which provides random-access to the source (by seeking to the given event ID or number)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "source = EventSource(get_dataset_path(\"gamma_prod5.simtel.zst\"), max_events=5)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Explore the contents of an event\n", - "\n", - "note that the R0 level is the raw data that comes out of a camera, and also the lowest level of monte-carlo data. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# so we can advance through events one-by-one\n", - "event_iterator = iter(source)\n", - "\n", - "event = next(event_iterator)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "the event is just a class with a bunch of data items in it. You can see a more compact represntation via:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "event.r0" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "printing the event structure, will currently print the value all items under it (so you get a lot of output if you print a high-level container):" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(event.simulation.shower)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(event.r0.tel.keys())" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "note that the event has 3 telescopes in it: Let's try the next one:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "event = next(event_iterator)\n", - "print(event.r0.tel.keys())" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "now, we have a larger event with many telescopes... Let's look at one of them:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "teldata = event.r0.tel[26]\n", - "print(teldata)\n", - "teldata" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that some values are unit quantities (`astropy.units.Quantity`) or angular quantities (`astropy.coordinates.Angle`), and you can easily maniuplate them:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "event.simulation.shower.energy" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "event.simulation.shower.energy.to(\"GeV\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "event.simulation.shower.energy.to(\"J\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "event.simulation.shower.alt" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(\"Altitude in degrees:\", event.simulation.shower.alt.deg)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Look for signal pixels in a camera\n", - "again, `event.r0.tel[x]` contains a data structure for the telescope data, with some fields like `waveform`.\n", - "\n", - "Let's make a 2D plot of the sample data (sample vs pixel), so we can see if we see which pixels contain Cherenkov light signals:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.pcolormesh(teldata.waveform[0]) # note the [0] is for channel 0\n", - "plt.colorbar()\n", - "plt.xlabel(\"sample number\")\n", - "plt.ylabel(\"Pixel_id\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's zoom in to see if we can identify the pixels that have the Cherenkov signal in them" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.pcolormesh(teldata.waveform[0])\n", - "plt.colorbar()\n", - "plt.ylim(700, 750)\n", - "plt.xlabel(\"sample number\")\n", - "plt.ylabel(\"pixel_id\")\n", - "print(\"waveform[0] is an array of shape (N_pix,N_slice) =\", teldata.waveform[0].shape)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now we can really see that some pixels have a signal in them!\n", - "\n", - "Lets look at a 1D plot of pixel 270 in channel 0 and see the signal:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "trace = teldata.waveform[0][719]\n", - "plt.plot(trace, drawstyle=\"steps\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Great! It looks like a *standard Cherenkov signal*!\n", - "\n", - "Let's take a look at several traces to see if the peaks area aligned:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for pix_id in range(718, 723):\n", - " plt.plot(\n", - " teldata.waveform[0][pix_id], label=\"pix {}\".format(pix_id), drawstyle=\"steps\"\n", - " )\n", - "plt.legend()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Look at the time trace from a Camera Pixel\n", - "\n", - "`ctapipe.calib.camera` includes classes for doing automatic trace integration with many methods, but before using that, let's just try to do something simple!\n", - "\n", - "Let's define the integration windows first:\n", - "By eye, they seem to be reaonsable from sample 8 to 13 for signal, and 20 to 29 for pedestal (which we define as the sum of all noise: NSB + electronic)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for pix_id in range(718, 723):\n", - " plt.plot(teldata.waveform[0][pix_id], \"+-\")\n", - "plt.fill_betweenx([0, 1600], 19, 24, color=\"red\", alpha=0.3, label=\"Ped window\")\n", - "plt.fill_betweenx([0, 1600], 5, 9, color=\"green\", alpha=0.3, label=\"Signal window\")\n", - "plt.legend()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Do a very simplisitic trace analysis \n", - "Now, let's for example calculate a signal and background in a the fixed windows we defined for this single event. Note we are ignoring the fact that cameras have 2 gains, and just using a single gain (channel 0, which is the high-gain channel):" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "data = teldata.waveform[0]\n", - "peds = data[:, 19:24].mean(axis=1) # mean of samples 20 to 29 for all pixels\n", - "sums = data[:, 5:9].sum(axis=1) / (13 - 8) # simple sum integration" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "phist = plt.hist(peds, bins=50, range=[0, 150])\n", - "plt.title(\"Pedestal Distribution of all pixels for a single event\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "let's now take a look at the pedestal-subtracted sums and a pedestal-subtracted signal:\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.plot(sums - peds)\n", - "plt.xlabel(\"pixel id\")\n", - "plt.ylabel(\"Pedestal-subtracted Signal\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now, we can clearly see that the signal is centered at 0 where there is no Cherenkov light, and we can also clearly see the shower around pixel 250." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# we can also subtract the pedestals from the traces themselves, which would be needed to compare peaks properly\n", - "for ii in range(270, 280):\n", - " plt.plot(data[ii] - peds[ii], drawstyle=\"steps\", label=\"pix{}\".format(ii))\n", - "plt.legend()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Camera Displays\n", - "\n", - "It's of course much easier to see the signal if we plot it in 2D with correct pixel positions! \n", - "\n", - ">note: the instrument data model is not fully implemented, so there is not a good way to load all the camera information (right now it is hacked into the `inst` sub-container that is read from the Monte-Carlo file)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "camgeom = source.subarray.tel[24].camera.geometry" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "title = \"CT24, run {} event {} ped-sub\".format(event.index.obs_id, event.index.event_id)\n", - "disp = CameraDisplay(camgeom, title=title)\n", - "disp.image = sums - peds\n", - "disp.cmap = plt.cm.RdBu_r\n", - "disp.add_colorbar()\n", - "disp.set_limits_percent(95) # autoscale" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "It looks like a nice signal! We have plotted our pedestal-subtracted trace integral, and see the shower clearly!\n", - "\n", - "Let's look at all telescopes:\n", - "\n", - "> note we plot here the raw signal, since we have not calculated the pedestals for each)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "for tel in event.r0.tel.keys():\n", - " plt.figure()\n", - " camgeom = source.subarray.tel[tel].camera.geometry\n", - " title = \"CT{}, run {} event {}\".format(\n", - " tel, event.index.obs_id, event.index.event_id\n", - " )\n", - " disp = CameraDisplay(camgeom, title=title)\n", - " disp.image = event.r0.tel[tel].waveform[0].sum(axis=1)\n", - " disp.cmap = plt.cm.RdBu_r\n", - " disp.add_colorbar()\n", - " disp.set_limits_percent(95)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## some signal processing...\n", - "\n", - "Let's try to detect the peak using the scipy.signal package:\n", - "https://docs.scipy.org/doc/scipy/reference/signal.html" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from scipy import signal\n", - "import numpy as np" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "pix_ids = np.arange(len(data))\n", - "has_signal = sums > 300\n", - "\n", - "widths = np.array(\n", - " [\n", - " 8,\n", - " ]\n", - ") # peak widths to search for (let's fix it at 8 samples, about the width of the peak)\n", - "peaks = [signal.find_peaks_cwt(trace, widths) for trace in data[has_signal]]\n", - "\n", - "for p, s in zip(pix_ids[has_signal], peaks):\n", - " print(\"pix{} has peaks at sample {}\".format(p, s))\n", - " plt.plot(data[p], drawstyle=\"steps-mid\")\n", - " plt.scatter(np.array(s), data[p, s])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "clearly the signal needs to be filtered first, or an appropriate wavelet used, but the idea is nice" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.8.13" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/docs/user-guide/tutorials/theta_square.ipynb b/docs/user-guide/tutorials/theta_square.ipynb deleted file mode 100644 index 6e2c7652f33..00000000000 --- a/docs/user-guide/tutorials/theta_square.ipynb +++ /dev/null @@ -1,212 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Make a theta-square plot\n", - "\n", - "This is a basic example to analyze some events and make a $\\Theta^2$ plot" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "ExecuteTime": { - "end_time": "2018-06-15T12:49:35.515499Z", - "start_time": "2018-06-15T12:49:34.968051Z" - } - }, - "outputs": [], - "source": [ - "%matplotlib inline" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "ExecuteTime": { - "end_time": "2018-06-15T12:49:37.807612Z", - "start_time": "2018-06-15T12:49:35.520552Z" - } - }, - "outputs": [], - "source": [ - "from astropy import units as u\n", - "from astropy.coordinates.angle_utilities import angular_separation\n", - "from astropy.coordinates import SkyCoord, AltAz\n", - "\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "\n", - "from ctapipe.io import EventSource\n", - "\n", - "from ctapipe.calib import CameraCalibrator\n", - "from ctapipe.image import ImageProcessor\n", - "from ctapipe.reco import ShowerProcessor\n", - "\n", - "from tqdm.auto import tqdm" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "ExecuteTime": { - "end_time": "2018-06-15T12:49:37.887391Z", - "start_time": "2018-06-15T12:49:37.818824Z" - } - }, - "source": [ - "Get source events in MC dataset. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "ExecuteTime": { - "end_time": "2018-06-15T12:49:37.887391Z", - "start_time": "2018-06-15T12:49:37.818824Z" - } - }, - "outputs": [], - "source": [ - "source = EventSource(\n", - " \"dataset://gamma_prod5.simtel.zst\",\n", - " # allowed_tels={1, 2, 3, 4},\n", - ")\n", - "\n", - "subarray = source.subarray\n", - "\n", - "calib = CameraCalibrator(subarray=subarray)\n", - "image_processor = ImageProcessor(subarray=subarray)\n", - "shower_processor = ShowerProcessor(subarray=subarray)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "ExecuteTime": { - "end_time": "2018-06-15T12:49:47.500199Z", - "start_time": "2018-06-15T12:49:37.893169Z" - } - }, - "outputs": [], - "source": [ - "off_angles = []\n", - "\n", - "for event in tqdm(source):\n", - "\n", - " # calibrating the event\n", - " calib(event)\n", - " image_processor(event)\n", - " shower_processor(event)\n", - "\n", - " reco_result = event.dl2.stereo.geometry[\"HillasReconstructor\"]\n", - "\n", - " # get angular offset between reconstructed shower direction and MC\n", - " # generated shower direction\n", - " true_shower = event.simulation.shower\n", - " off_angle = angular_separation(\n", - " true_shower.az, true_shower.alt, reco_result.az, reco_result.alt\n", - " )\n", - "\n", - " # Appending all estimated off angles\n", - " off_angles.append(off_angle.to(u.deg).value)" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "ExecuteTime": { - "end_time": "2018-06-15T12:49:47.507369Z", - "start_time": "2018-06-15T12:49:47.502642Z" - } - }, - "source": [ - "calculate theta square for angles which are not nan" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "ExecuteTime": { - "end_time": "2018-06-15T12:49:47.507369Z", - "start_time": "2018-06-15T12:49:47.502642Z" - } - }, - "outputs": [], - "source": [ - "off_angles = np.array(off_angles)\n", - "thetasquare = off_angles[np.isfinite(off_angles)] ** 2" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "ExecuteTime": { - "end_time": "2018-06-15T12:49:48.264122Z", - "start_time": "2018-06-15T12:49:47.511172Z" - } - }, - "source": [ - "## Plot the results" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "ExecuteTime": { - "end_time": "2018-06-15T12:49:48.264122Z", - "start_time": "2018-06-15T12:49:47.511172Z" - } - }, - "outputs": [], - "source": [ - "plt.hist(thetasquare, bins=10, range=[0, 0.4])\n", - "plt.xlabel(r\"$\\theta^2$ (deg)\")\n", - "plt.ylabel(\"# of events\")\n", - "plt.show()" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.8.13" - }, - "toc": { - "nav_menu": { - "height": "13px", - "width": "253px" - }, - "number_sections": false, - "sideBar": true, - "skip_h1_title": false, - "toc_cell": false, - "toc_position": {}, - "toc_section_display": "block", - "toc_window_display": false - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/environment.yml b/environment.yml index e0951f7fec1..7f0f0dadeed 100644 --- a/environment.yml +++ b/environment.yml @@ -22,6 +22,7 @@ dependencies: - numpy>=1.22 - numpydoc - pandas + - pypandoc - pre-commit - psutil - pytables @@ -34,6 +35,7 @@ dependencies: - scipy - setuptools - sphinx + - sphinx-gallery - sphinx-automodapi - pydata-sphinx-theme - sphinx-design diff --git a/examples/README.rst b/examples/README.rst new file mode 100644 index 00000000000..cdc806953ae --- /dev/null +++ b/examples/README.rst @@ -0,0 +1,6 @@ +.. _tutorials_and_examples_gallery: + +Tutorials and Examples +====================== + +This gallery provides an overview of different ctapipe modules and how to use them. diff --git a/examples/algorithms/README.rst b/examples/algorithms/README.rst new file mode 100644 index 00000000000..a53d619994e --- /dev/null +++ b/examples/algorithms/README.rst @@ -0,0 +1,4 @@ +.. _algorithms-examples-gallery: + +Algorithms +---------- diff --git a/examples/algorithms/convert_images_to_2d.py b/examples/algorithms/convert_images_to_2d.py new file mode 100644 index 00000000000..9903df780f1 --- /dev/null +++ b/examples/algorithms/convert_images_to_2d.py @@ -0,0 +1,96 @@ +""" +Convert camera images to pixels on a s square grid +================================================== + +""" + +import astropy.units as u +import matplotlib.pyplot as plt + +from ctapipe.image.toymodel import Gaussian +from ctapipe.instrument import SubarrayDescription +from ctapipe.visualization import CameraDisplay + +###################################################################### +# get the subarray from an example file +subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst") + + +###################################################################### +# Geometries with square pixels +# ----------------------------- +# +# Define a camera geometry and generate a dummy image: +# + +geom = subarray.tel[40].camera.geometry +model = Gaussian( + x=0.05 * u.m, + y=0.05 * u.m, + width=0.01 * u.m, + length=0.05 * u.m, + psi="30d", +) +_, image, _ = model.generate_image(geom, intensity=500, nsb_level_pe=3) + +###################################################################### +CameraDisplay(geom, image) + + +###################################################################### +# The ``CameraGeometry`` has functions to convert the 1d image arrays to +# 2d arrays and back to the 1d array: +# + +image_square = geom.image_to_cartesian_representation(image) + +###################################################################### +plt.imshow(image_square) + +###################################################################### +image_1d = geom.image_from_cartesian_representation(image_square) + +###################################################################### +CameraDisplay(geom, image_1d) + + +###################################################################### +# Geometries with hexagonal pixels +# -------------------------------- +# +# Define a camera geometry and generate a dummy image: +# + +geom = subarray.tel[1].camera.geometry +model = Gaussian( + x=0.5 * u.m, + y=0.5 * u.m, + width=0.1 * u.m, + length=0.2 * u.m, + psi="30d", +) +_, image, _ = model.generate_image(geom, intensity=5000) + +###################################################################### +CameraDisplay(geom, image) + +###################################################################### +image_square = geom.image_to_cartesian_representation(image) + + +###################################################################### +# Conversion into square geometry +# ------------------------------- +# +# Since the resulting array has square pixels, the pixel grid has to be +# rotated and distorted. This is reversible (The +# ``image_from_cartesian_representation`` method takes care of this): +# + +plt.imshow(image_square) + +###################################################################### +image_1d = geom.image_from_cartesian_representation(image_square) + +###################################################################### +disp = CameraDisplay(geom, image_1d) diff --git a/examples/algorithms/dilate_image.py b/examples/algorithms/dilate_image.py new file mode 100644 index 00000000000..a5757be43b6 --- /dev/null +++ b/examples/algorithms/dilate_image.py @@ -0,0 +1,67 @@ +""" +Basic Image Cleaning and Dilation +================================= + +Here we create an example shower image, do a tail-cuts +(picture/boundary) cleaning, and then dilate the resulting cleaning mask +by several neighbor pixels + +""" + +import astropy.units as u + +# %matplotlib inline +from matplotlib import pyplot as plt + +from ctapipe.image import dilate, tailcuts_clean, toymodel +from ctapipe.instrument import SubarrayDescription +from ctapipe.visualization import CameraDisplay + +# Load a camera from an example file +subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst") +geom = subarray.tel[100].camera.geometry + +# Create a fake camera image to display: +model = toymodel.Gaussian( + x=0.2 * u.m, + y=0.0 * u.m, + width=0.05 * u.m, + length=0.15 * u.m, + psi="35d", +) + +image, sig, bg = model.generate_image(geom, intensity=1500, nsb_level_pe=5) + + +###################################################################### +# Apply the image cleaning: +# + +cleanmask = tailcuts_clean(geom, image, picture_thresh=10, boundary_thresh=5) +clean = image.copy() +clean[~cleanmask] = 0.0 + +disp = CameraDisplay(geom, image=image) +disp.highlight_pixels(cleanmask, color="red") + + +###################################################################### +# Now dialte the mask a few times: +# + + +def show_dilate(mask, times=1): + m = mask.copy() + for ii in range(times): + m = dilate(geom, m) + CameraDisplay( + geom, image=(m.astype(int) + mask.astype(int)), title="dilate{}".format(times) + ) + + +###################################################################### +plt.figure(figsize=(18, 3)) + +for ii in range(0, 6): + plt.subplot(1, 7, ii + 1) + show_dilate(cleanmask.copy(), times=ii) diff --git a/examples/algorithms/nd_interpolation.py b/examples/algorithms/nd_interpolation.py new file mode 100644 index 00000000000..21a03057e1a --- /dev/null +++ b/examples/algorithms/nd_interpolation.py @@ -0,0 +1,144 @@ +""" +Use N-dimensional Histogram functionality and Interpolation +=========================================================== + +- could be used for example to read and interpolate an lookup table or + IRF. +- In this example, we load a sample energy reconstruction lookup-table + from a FITS file +- In this case it is only in 2D cube (to keep the file size small): + ``SIZE`` vs ``IMPACT-DISTANCE``, however the same method will work + for any dimensionality + +""" + +import matplotlib.pylab as plt +import numpy as np +from astropy.io import fits +from scipy.interpolate import RegularGridInterpolator + +from ctapipe.utils import Histogram +from ctapipe.utils.datasets import get_dataset_path + +# %matplotlib inline + + +###################################################################### +# load an example datacube +# ~~~~~~~~~~~~~~~~~~~~~~~~ +# +# (an energy table generated for a small subset of HESS simulations) to +# use as a lookup table. Here we will use the ``Histogram`` class, which +# automatically loads both the data cube and creates arrays for the +# coordinates of each axis. +# + +testfile = get_dataset_path("hess_ct001_energylut.fits.gz") +energy_hdu = fits.open(testfile)["MEAN"] +energy_table = Histogram.from_fits(energy_hdu) +print(energy_table) + + +###################################################################### +# construct an interpolator that we can use to get values at any point: +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Here we will use a ``RegularGridInterpolator``, since it is the most +# appropriate for this type of data, but others are available (see the +# SciPy documentation) +# + +centers = [energy_table.bin_centers(ii) for ii in range(energy_table.ndims)] +energy_lookup = RegularGridInterpolator( + centers, energy_table.hist, bounds_error=False, fill_value=-100 +) + + +###################################################################### +# ``energy_lookup`` is now just a continuous function of ``log(SIZE)``, +# ``DISTANCE`` in m. +# +# Now plot some curves from the interpolator. +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Note that the LUT we used is does not have very high statistics, so the +# interpolation starts to be affected by noise at the high end. In a real +# case, we would want to use a table that has been sanitized (smoothed and +# extrapolated) +# + +lsize = np.linspace(1.5, 5.0, 100) +dists = np.linspace(50, 100, 5) + +plt.figure() +plt.title("Variation of energy with size and impact distance") +plt.xlabel("SIZE (P.E.)") +plt.ylabel("ENERGY (TeV)") + +for dist in dists: + plt.plot( + 10**lsize, + 10 ** energy_lookup((lsize, dist)), + "+-", + label="DIST={:.1f} m".format(dist), + ) + +plt.legend(loc="best") + + +###################################################################### +# Using the interpolator, reinterpolate the lookup table onto an +# :math:`N \times N` grid (regardless of its original dimensions): +# + +N = 300 +xmin, xmax = energy_table.bin_centers(0)[0], energy_table.bin_centers(0)[-1] +ymin, ymax = energy_table.bin_centers(1)[0], energy_table.bin_centers(1)[-1] +xx, yy = np.linspace(xmin, xmax, N), np.linspace(ymin, ymax, N) +X, Y = np.meshgrid(xx, yy) +points = list(zip(X.ravel(), Y.ravel())) +E = energy_lookup(points).reshape((N, N)) + + +###################################################################### +# Now, let’s plot the original table and the new one (E). The color bar +# shows :math:`\log_{10}(E)` in TeV +# + +plt.figure(figsize=(12, 5)) +plt.nipy_spectral() + +# the uninterpolated table +plt.subplot(1, 2, 1) +plt.xlim(1.5, 5) +plt.ylim(0, 500) +plt.xlabel("log10(SIZE)") +plt.ylabel("Impact Dist (m)") +plt.pcolormesh( + energy_table.bin_centers(0), energy_table.bin_centers(1), energy_table.hist.T +) +plt.title("Raw table, uninterpolated {0}".format(energy_table.hist.T.shape)) +cb = plt.colorbar() +cb.set_label(r"$\log_{10}(E/\mathrm{TeV})$") + +# the interpolated table +plt.subplot(1, 2, 2) +plt.pcolormesh(np.linspace(xmin, xmax, N), np.linspace(ymin, ymax, N), E) +plt.xlim(1.5, 5) +plt.ylim(0, 500) +plt.xlabel("log10(SIZE)") +plt.ylabel("Impact Dist(m)") +plt.title("Interpolated to a ({0}, {0}) grid".format(N)) +cb = plt.colorbar() +cb.set_label(r"$\log_{10}(E/\mathrm{TeV})$") + +plt.tight_layout() +plt.show() + + +###################################################################### +# In the high-stats central region, we get a nice smooth interpolation +# function. Of course we can see that there are a few more steps to take +# before using this table: \* +# - need to deal with cases where the table had low stats near the edges (smooth or extrapolate, or set bounds) +# - may need to smooth the table even where there are sufficient stats, to avoid wiggles diff --git a/examples/core/InstrumentDescription.py b/examples/core/InstrumentDescription.py new file mode 100644 index 00000000000..6d26ab7f9b3 --- /dev/null +++ b/examples/core/InstrumentDescription.py @@ -0,0 +1,174 @@ +""" +Working with Instrumental Descriptions +====================================== + +the instrumental description is loaded by the event source, and consists +of a hierarchy of classes in the ctapipe.instrument module, the base of +which is the ``SubarrayDescription`` + +""" + +from astropy.coordinates import SkyCoord + +from ctapipe.io import EventSource +from ctapipe.utils.datasets import get_dataset_path +from ctapipe.visualization import CameraDisplay + +filename = get_dataset_path("gamma_prod5.simtel.zst") + +with EventSource(filename, max_events=1) as source: + subarray = source.subarray + + +###################################################################### +# the SubarrayDescription: +# ------------------------ +# + +subarray.info() + +###################################################################### +subarray.to_table() + + +###################################################################### +# You can also get a table of just the ``OpticsDescriptions`` +# (``CameraGeometry`` is more complex and can’t be stored on a single +# table row, so each one can be converted to a table separately) +# + +subarray.to_table(kind="optics") + + +###################################################################### +# Make a sub-array with only SC-type telescopes: +# + +sc_tels = [tel_id for tel_id, tel in subarray.tel.items() if tel.optics.n_mirrors == 2] +newsub = subarray.select_subarray(sc_tels, name="SCTels") +newsub.info() + + +###################################################################### +# can also do this by using ``Table.group_by`` +# + + +###################################################################### +# Explore some of the details of the telescopes +# --------------------------------------------- +# + +tel = subarray.tel[1] +tel + +###################################################################### +tel.optics.mirror_area + +###################################################################### +tel.optics.n_mirror_tiles + +###################################################################### +tel.optics.equivalent_focal_length + +###################################################################### +tel.camera + +###################################################################### +tel.camera.geometry.pix_x + +###################################################################### +# %matplotlib inline + +CameraDisplay(tel.camera.geometry) + +###################################################################### +CameraDisplay(subarray.tel[98].camera.geometry) + + +###################################################################### +# Plot the subarray +# ----------------- +# +# We’ll make a subarray by telescope type and plot each separately, so +# they appear in different colors. We also calculate the radius using the +# mirror area (and exagerate it a bit). +# +# This is just for debugging and info, for any “real” use, a +# ``visualization.ArrayDisplay`` should be used +# + +subarray.peek() + +###################################################################### +subarray.footprint + + +###################################################################### +# Get info about the subarray in general +# -------------------------------------- +# + +subarray.telescope_types + +###################################################################### +subarray.camera_types + +###################################################################### +subarray.optics_types + + +###################################################################### +center = SkyCoord("10.0 m", "2.0 m", "0.0 m", frame="groundframe") +coords = subarray.tel_coords # a flat list of coordinates by tel_index +coords.separation(center) + + +###################################################################### +# Telescope IDs vs Indices +# ------------------------ +# +# Note that ``subarray.tel`` is a dict mapped by ``tel_id`` (the +# indentifying number of a telescope). It is possible to have telescope +# IDs that do not start at 0, are not contiguouous (e.g. if a subarray is +# selected). Some functions and properties like ``tel_coords`` are numpy +# arrays (not dicts) so they are not mapped to the telescope ID, but +# rather the *index* within this SubarrayDescription. To convert between +# the two concepts you can do: +# + +subarray.tel_ids_to_indices([1, 5, 23]) + + +###################################################################### +# or you can get the indexing array directly in numpy or dict form: +# + +subarray.tel_index_array + +###################################################################### +subarray.tel_index_array[[1, 5, 23]] + +###################################################################### +subarray.tel_indices[ + 1 +] # this is a dict of tel_id -> tel_index, so we can only do one at once + +ids = subarray.get_tel_ids_for_type(subarray.telescope_types[0]) +ids + +###################################################################### +idx = subarray.tel_ids_to_indices(ids) +idx + +###################################################################### +subarray.tel_coords[idx] + + +###################################################################### +# so, with that method you can quickly get many telescope positions at +# once (the alternative is to use the dict ``positions`` which maps +# ``tel_id`` to a position on the ground +# + +subarray.positions[1] diff --git a/examples/core/README.rst b/examples/core/README.rst new file mode 100644 index 00000000000..cfd78a64e9e --- /dev/null +++ b/examples/core/README.rst @@ -0,0 +1,4 @@ +.. _core-examples-gallery: + +Core Functionality +------------------ diff --git a/examples/core/command_line_tools.py b/examples/core/command_line_tools.py new file mode 100644 index 00000000000..195a8e880a1 --- /dev/null +++ b/examples/core/command_line_tools.py @@ -0,0 +1,359 @@ +""" +Creating command-line Tools +=========================== + +""" + +from time import sleep + +from ctapipe.core import Component, TelescopeComponent, Tool +from ctapipe.core.traits import ( + FloatTelescopeParameter, + Integer, + Path, + TraitError, + observe, +) +from ctapipe.instrument import SubarrayDescription +from ctapipe.utils import get_dataset_path + +###################################################################### +GAMMA_FILE = get_dataset_path("gamma_prod5.simtel.zst") + + +###################################################################### +# see https://github.com/ipython/traitlets/blob/master/examples/myapp.py +# + + +###################################################################### +# Setup: +# ------ +# +# Create a few ``Component``\ s that we will use later in a ``Tool``: +# + + +class MyComponent(Component): + """A Component that does stuff""" + + value = Integer(default_value=-1, help="Value to use").tag(config=True) + + def do_thing(self): + self.log.debug("Did thing") + + +# in order to have 2 of the same components at once +class SecondaryMyComponent(MyComponent): + """A second component""" + + +class AdvancedComponent(Component): + """An advanced technique""" + + value1 = Integer(default_value=-1, help="Value to use").tag(config=True) + infile = Path( + help="input file name", + exists=None, # set to True to require existing, False for requiring non-existing + directory_ok=False, + ).tag(config=True) + outfile = Path(help="output file name", exists=False, directory_ok=False).tag( + config=True + ) + + def __init__(self, config=None, parent=None, **kwargs): + super().__init__(config=config, parent=parent, **kwargs) + # components can have sub components, but these must have + # then parent=self as argument and be assigned as member + # so the full config can be received later + self.subcompent = MyComponent(parent=self) + + @observe("outfile") + def on_outfile_changed(self, change): + self.log.warning("Outfile was changed to '{}'".format(change)) + + +class TelescopeWiseComponent(TelescopeComponent): + """a component that contains parameters that are per-telescope configurable""" + + param = FloatTelescopeParameter( + help="Something configurable with telescope patterns", default_value=5.0 + ).tag(config=True) + + +###################################################################### +MyComponent() + +###################################################################### +AdvancedComponent(infile="test.foo", outfile="out.foo") + + +###################################################################### +# ``TelescopeComponents`` need to have a subarray given to them in order +# to work (since they need one to turn a ``TelescopeParameter`` into a +# concrete list of values for each telescope. Here we will give a dummy +# one: +# + + +subarray = SubarrayDescription.read(GAMMA_FILE) +subarray.info() + +###################################################################### +TelescopeWiseComponent(subarray=subarray) + + +###################################################################### +# This TelescopeParameters can then be set using a list of patterns like: +# +# +# component.param = [("type", "LST*",3.0),("type", "MST*", 2.0),(id, 25, 4.0)] +# +# These get translated into per-telescope-id values once the subarray is +# registered. After that one acccess the per-telescope id values via: +# +# +# component.param.tel[tel_id] +# + + +###################################################################### +# Now create an executable Tool that contains the Components +# ---------------------------------------------------------- +# +# Note that all the components we wish to be configured via the tool must +# be added to the ``classes`` attribute. +# + + +class MyTool(Tool): + name = "mytool" + description = "do some things and stuff" + aliases = dict( + infile="AdvancedComponent.infile", + outfile="AdvancedComponent.outfile", + iterations="MyTool.iterations", + ) + + # Which classes are registered for configuration + classes = [ + MyComponent, + AdvancedComponent, + SecondaryMyComponent, + TelescopeWiseComponent, + ] + + # local configuration parameters + iterations = Integer(5, help="Number of times to run", allow_none=False).tag( + config=True + ) + + def setup(self): + self.comp = MyComponent(parent=self) + self.comp2 = SecondaryMyComponent(parent=self) + self.comp3 = TelescopeWiseComponent(parent=self, subarray=subarray) + self.advanced = AdvancedComponent(parent=self) + + def start(self): + self.log.info("Performing {} iterations...".format(self.iterations)) + for ii in range(self.iterations): + self.log.info("ITERATION {}".format(ii)) + self.comp.do_thing() + self.comp2.do_thing() + sleep(0.1) + + def finish(self): + self.log.warning("Shutting down.") + + +###################################################################### +# Get Help info +# ------------- +# +# The following allows you to print the help info within a Jupyter +# notebook, but this same inforamtion would be displayed if the user +# types: +# +# :: +# +# mytool --help +# + +tool = MyTool() +tool + +###################################################################### +tool.print_help() + + +###################################################################### +# The following is equivalant to the user typing ``mytool --help-all`` +# + +tool.print_help(classes=True) + + +###################################################################### +# Run the tool +# ------------ +# +# here we pass in argv since it is a Notebook, but if argv is not +# specified it’s read from ``sys.argv``, so the following is the same as +# running: +# +# +# mytool --log_level=INFO --infile gamma_test.simtel.gz --iterations=3 +# +# As Tools are intended to be exectutables, they are raising +# ``SystemExit`` on exit. Here, we use them to demonstrate how it would +# work, so we catch the ``SystemExit``. +# + +try: + tool.run(argv=["--infile", str(GAMMA_FILE), "--outfile", "out.csv"]) +except SystemExit as e: + assert e.code == 0, f"Tool returned with error status {e}" + +tool.log_format = "%(asctime)s : %(levelname)s [%(name)s %(funcName)s] %(message)s" + + +###################################################################### +try: + tool.run( + argv=[ + "--log-level", + "INFO", + "--infile", + str(GAMMA_FILE), + "--outfile", + "out.csv", + "--iterations", + "3", + ] + ) +except SystemExit as e: + assert e.code == 0, f"Tool returned with error status {e}" + + +###################################################################### +# here we change the log-level to DEBUG: +# + +try: + tool.run( + argv=[ + "--log-level", + "DEBUG", + "--infile", + str(GAMMA_FILE), + "--outfile", + "out.csv", + ] + ) +except SystemExit as e: + assert e.code == 0, f"Tool returned with error status {e}" + + +###################################################################### +# you can also set parameters directly in the class, rather than using the +# argument/configfile parser. This is useful if you are calling the Tool +# from a script rather than the command-line +# + +tool.iterations = 1 +tool.log_level = 0 + +try: + tool.run(["--infile", str(GAMMA_FILE), "--outfile", "out.csv"]) +except SystemExit as e: + assert e.code == 0, f"Tool returned with error status {e}" + + +###################################################################### +# see what happens when a value is set that is not of the correct type: +# + +try: + tool.iterations = "badval" +except TraitError as E: + print("bad value:", E) +except SystemExit as e: + assert e.code == 0, f"Tool returned with error status {e}" + + +###################################################################### +# Example of what happens when you change a parameter that is being +# “observed” in a class. It’s handler is called: +# + +tool.advanced.outfile = "Another.txt" + + +###################################################################### +# we see that the handler for ``outfile`` was called, and it receive a +# change dict that shows the old and new values. +# + + +###################################################################### +# create a tool using a config file: +# + +tool2 = MyTool() + +###################################################################### +try: + tool2.run(argv=["--config", "config.json"]) +except SystemExit as e: + assert e.code == 0, f"Tool returned with error status {e}" + +###################################################################### +print(tool2.advanced.infile) + +###################################################################### +print(tool2.config) + +###################################################################### +tool2.is_setup + +###################################################################### +tool3 = MyTool() + +###################################################################### +tool3.is_setup + +###################################################################### +tool3.initialize(argv=[]) + +###################################################################### +tool3.is_setup + +###################################################################### +tool3 + +###################################################################### +tool.setup() +tool + +###################################################################### +tool.comp2 + + +###################################################################### +# Getting the configuration of an instance +# ---------------------------------------- +# + +tool.get_current_config() + +###################################################################### +tool.iterations = 12 +tool.get_current_config() + + +###################################################################### +# Writing a Sample Config File +# ---------------------------- +# + +print(tool.generate_config_file()) diff --git a/docs/user-guide/examples/Tools.json b/examples/core/config.json similarity index 100% rename from docs/user-guide/examples/Tools.json rename to examples/core/config.json diff --git a/examples/core/containers.py b/examples/core/containers.py new file mode 100644 index 00000000000..82242604b55 --- /dev/null +++ b/examples/core/containers.py @@ -0,0 +1,212 @@ +""" +Using Container classes +======================= + +``ctapipe.core.Container`` is the base class for all event-wise data +classes in ctapipe. It works like a object-relational mapper, in that it +defines a set of ``Fields`` along with their metadata (description, +unit, default), which can be later translated automatially into an +output table using a ``ctapipe.io.TableWriter``. + +""" + +from functools import partial + +import numpy as np +from astropy import units as u + +from ctapipe.containers import SimulatedShowerContainer +from ctapipe.core import Container, Field, Map + +###################################################################### +# Let’s define a few example containers with some dummy fields in them: +# + + +class SubContainer(Container): + junk = Field(-1, "Some junk") + value = Field(0.0, "some value", unit=u.deg) + + +class TelContainer(Container): + # defaults should match the other requirements, e.g. the defaults + # should have the correct unit. It most often also makes sense to use + # an invalid value marker like nan for floats or -1 for positive integers + # as default + tel_id = Field(-1, "telescope ID number") + + # For mutable structures like lists, arrays or containers, use a `default_factory` function or class + # not an instance to assure each container gets a fresh instance and there is no hidden + # shared state between containers. + image = Field(default_factory=lambda: np.zeros(10), description="camera pixel data") + + +class EventContainer(Container): + event_id = Field(-1, "event id number") + + tels_with_data = Field( + default_factory=list, description="list of telescopes with data" + ) + sub = Field( + default_factory=SubContainer, description="stuff" + ) # a sub-container in the hierarchy + + # A Map is like a defaultdictionary with a specific container type as default. + # This can be used to e.g. store a container per telescope + # we use partial here to automatically get a function that creates a map with the correct container type + # as default + tel = Field(default_factory=partial(Map, TelContainer), description="telescopes") + + +###################################################################### +# Basic features +# -------------- +# + +ev = EventContainer() + + +###################################################################### +# Check that default values are automatically filled in +# + +print(ev.event_id) +print(ev.sub) +print(ev.tel) +print(ev.tel.keys()) + +# default dict access will create container: +print(ev.tel[1]) + + +###################################################################### +# print the dict representation +# + +print(ev) + + +###################################################################### +# We also get docstrings “for free” +# +help(EventContainer) + +###################################################################### +help(SubContainer) + +###################################################################### +# values can be set as normal for a class: +# + +ev.event_id = 100 +ev.event_id + +###################################################################### +ev.as_dict() # by default only shows the bare items, not sub-containers (See later) + +###################################################################### +ev.as_dict(recursive=True) + + +###################################################################### +# and we can add a few of these to the parent container inside the tel +# dict: +# + +ev.tel[10] = TelContainer() +ev.tel[5] = TelContainer() +ev.tel[42] = TelContainer() + +###################################################################### +# because we are using a default_factory to handle mutable defaults, the images are actually different: +ev.tel[42].image is ev.tel[32] + + +###################################################################### +# Be careful to use the ``default_factory`` mechanism for mutable fields, +# see this **negative** example: +# + + +class DangerousContainer(Container): + image = Field( + np.zeros(10), + description=( + "Attention!!!! Globally mutable shared state. Use default_factory instead" + ), + ) + + +c1 = DangerousContainer() +c2 = DangerousContainer() + +c1.image[5] = 9999 + +print(c1.image) +print(c2.image) +print(c1.image is c2.image) + +###################################################################### +ev.tel + + +###################################################################### +# Converion to dictionaries +# ------------------------- +# + +ev.as_dict() + +###################################################################### +ev.as_dict(recursive=True, flatten=False) + + +###################################################################### +# for serialization to a table, we can even flatten the output into a +# single set of columns +# + +ev.as_dict(recursive=True, flatten=True) + + +###################################################################### +# Setting and clearing values +# --------------------------- +# + +ev.tel[5].image[:] = 9 +print(ev) + +###################################################################### +ev.reset() +ev.as_dict(recursive=True) + + +###################################################################### +# look at a pre-defined Container +# ------------------------------- +# + + +help(SimulatedShowerContainer) + +###################################################################### +shower = SimulatedShowerContainer() +shower + + +###################################################################### +# Container prefixes +# ------------------ +# +# To store the same container in the same table in a file or give more +# information, containers support setting a custom prefix: +# + +c1 = SubContainer(junk=5, value=3, prefix="foo") +c2 = SubContainer(junk=10, value=9001, prefix="bar") + +# create a common dict with data from both containers: +d = c1.as_dict(add_prefix=True) +d.update(c2.as_dict(add_prefix=True)) +d diff --git a/examples/core/provenance.py b/examples/core/provenance.py new file mode 100644 index 00000000000..2a5ae103bfe --- /dev/null +++ b/examples/core/provenance.py @@ -0,0 +1,118 @@ +""" +Using the ctapipe Provenance service +==================================== + +The provenance functionality is used automatically when you use most of +ctapipe functionality (particularly ``ctapipe.core.Tool`` and functions +in ``ctapipe.io`` and ``ctapipe.utils``), so normally you don’t have to +work with it directly. It tracks both input and output files, as well as +details of the machine and software environment on which a Tool +executed. + +Here we show some very low-level functions of this system: + +""" + +from pprint import pprint + +from ctapipe.core import Provenance + +###################################################################### +# Activities +# ---------- +# +# The basis of Provenance is an *activity*, which is generally an +# executable or step in a script. Activities can be nested (e.g. with +# sub-activities), as shown below, but normally this is not required: +# + +p = Provenance() # note this is a singleton, so only ever one global provenence object +p.clear() +p.start_activity() +p.add_input_file("test.txt") + +p.start_activity("sub") +p.add_input_file("subinput.txt") +p.add_input_file("anothersubinput.txt") +p.add_output_file("suboutput.txt") +p.finish_activity("sub") + +p.start_activity("sub2") +p.add_input_file("sub2input.txt") +p.finish_activity("sub2") + +p.finish_activity() + +###################################################################### +p.finished_activity_names + + +###################################################################### +# Activities have associated input and output *entities* (files or other +# objects) +# + +[(x["activity_name"], x["input"]) for x in p.provenance] + + +###################################################################### +# Activities track when they were started and finished: +# + +[(x["activity_name"], x["duration_min"]) for x in p.provenance] + + +###################################################################### +# Full provenance +# --------------- +# +# The provence object is a list of activitites, and for each lots of +# details are collected: +# + +p.provenance[0] + + +###################################################################### +# This can be better represented in JSON: +# + +print(p.as_json(indent=2)) + + +###################################################################### +# Storing provenance info in output files +# --------------------------------------- +# +# - already this can be stored in something like an HDF5 file header, +# which allows hierarchies. +# - Try to flatted the data so it can be stored in a key=value header in +# a **FITS file** (using the FITS extended keyword convention to allow +# >8 character keywords), or as a table +# + + +def flatten_dict(y): + out = {} + + def flatten(x, name=""): + if type(x) is dict: + for a in x: + flatten(x[a], name + a + ".") + elif type(x) is list: + i = 0 + for a in x: + flatten(a, name + str(i) + ".") + i += 1 + else: + out[name[:-1]] = x + + flatten(y) + return out + + +###################################################################### +d = dict(activity=p.provenance) + +###################################################################### +pprint(flatten_dict(d)) diff --git a/examples/core/table_writer_reader.py b/examples/core/table_writer_reader.py new file mode 100644 index 00000000000..e578c264123 --- /dev/null +++ b/examples/core/table_writer_reader.py @@ -0,0 +1,289 @@ +""" +Writing Containers to a tabular format +====================================== + +The ``TableWriter``/``TableReader`` sub-classes allow you to write a +``ctapipe.core.Container`` class and its meta-data to an output table. +They treat the ``Field``s in the ``Container`` as columns in the +output, and automatically generate a schema. Here we will go through an +example of writing out data and reading it back with *Pandas*, +*PyTables*, and a ``ctapipe.io.TableReader``: + +In this example, we will use the ``HDF5TableWriter``, which writes to +HDF5 datasets using *PyTables*. Currently this is the only implemented +TableWriter. + +""" + + +###################################################################### +# Caveats to think about: \* vector columns in Containers *can* be +# written, but some lilbraries like Pandas can not read those (so you must +# use pytables or astropy to read outputs that have vector columns) +# \* units are stored in the table metadata, but some libraries like Pandas +# ignore them and all other metadata +# + + +###################################################################### +# Create some example Containers +# ------------------------------ +# + +import os + +import numpy as np +import pandas as pd +import tables +from astropy import units as u + +from ctapipe.core import Container, Field +from ctapipe.io import HDF5TableReader, HDF5TableWriter, read_table + + +###################################################################### +class VariousTypesContainer(Container): + + a_int = Field(int, "some int value") + a_float = Field(float, "some float value with a unit", unit=u.m) + a_bool = Field(bool, "some bool value") + a_np_int = Field(np.int64, "a numpy int") + a_np_float = Field(np.float64, "a numpy float") + a_np_bool = Field(np.bool_, "np.bool") + + +###################################################################### +# let’s also make a dummy stream (generator) that will create a series of +# these containers +# + + +def create_stream(n_event): + + data = VariousTypesContainer() + for i in range(n_event): + + data.a_int = int(i) + data.a_float = float(i) * u.cm # note unit conversion will happen + data.a_bool = (i % 2) == 0 + data.a_np_int = np.int64(i) + data.a_np_float = np.float64(i) + data.a_np_bool = np.bool_((i % 2) == 0) + + yield data + + +###################################################################### +for data in create_stream(2): + + for key, val in data.items(): + + print("{}: {}, type : {}".format(key, val, type(val))) + + +###################################################################### +# Writing the Data (and good practices) +# ------------------------------------- +# + + +###################################################################### +# Always use context managers with IO classes, as they will make sure the +# underlying resources are properly closed in case of errors: +# + +try: + with HDF5TableWriter("container.h5", group_name="data") as h5_table: + + for data in create_stream(10): + + h5_table.write("table", data) + 0 / 0 +except Exception as err: + print("FAILED:", err) +print("Done") + +h5_table.h5file.isopen == False + +###################################################################### +print(os.listdir()) + +###################################################################### +# Appending new Containers +# ------------------------ +# + + +###################################################################### +# To append some new containers we need to set the writing in append mode +# by using: ‘mode=a’. But let’s now first look at what happens if we +# don’t. +# + +for i in range(2): + + with HDF5TableWriter( + "container.h5", mode="w", group_name="data_{}".format(i) + ) as h5_table: + + for data in create_stream(10): + + h5_table.write("table", data) + + print(h5_table.h5file) + +###################################################################### +os.remove("container.h5") + +###################################################################### +# Ok so the writer destroyed the content of the file each time it opens +# the file. Now let’s try to append some data group to it! (using +# mode=‘a’) +# + +for i in range(2): + + with HDF5TableWriter( + "container.h5", mode="a", group_name="data_{}".format(i) + ) as h5_table: + + for data in create_stream(10): + + h5_table.write("table", data) + + print(h5_table.h5file) + + +###################################################################### +# So we can append some data groups. As long as the data group_name does +# not already exists. Let’s try to overwrite the data group : data_1 +# + +try: + with HDF5TableWriter("container.h5", mode="a", group_name="data_1") as h5_table: + for data in create_stream(10): + h5_table.write("table", data) +except Exception as err: + print("Failed as expected:", err) + + +###################################################################### +# Good ! I cannot overwrite my data. +# + +print(bool(h5_table.h5file.isopen)) + + +###################################################################### +# Reading the Data +# ---------------- +# + + +###################################################################### +# Reading the whole table at once: +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# For this, you have several choices. Since we used the HDF5TableWriter in +# this example, we have at least these options avilable: +# +# - Pandas +# - PyTables +# - Astropy Table +# +# For other TableWriter implementations, others may be possible (depending +# on format) +# + + +###################################################################### +# Reading using ``ctapipe.io.read_table`` +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +# This is the preferred method, it returns an astropy ``Table`` and +# supports keeping track of units, metadata and transformations. +# + + +table = read_table("container.h5", "/data_0/table") +table[:5] + +###################################################################### +table.meta + + +###################################################################### +# Reading with Pandas: +# ^^^^^^^^^^^^^^^^^^^^ +# +# Pandas is a convenient way to read the output. **HOWEVER BE WARNED** +# that so far Pandas does not support reading the table *meta-data* or +# *units* for colums, so that information is lost! +# + + +data = pd.read_hdf("container.h5", key="/data_0/table") +data.head() + + +###################################################################### +# Reading with PyTables +# ^^^^^^^^^^^^^^^^^^^^^ +# + + +h5 = tables.open_file("container.h5") +table = h5.root["data_0"]["table"] +table + + +###################################################################### +# note that here we can still access the metadata +# + +table.attrs + + +###################################################################### +# Reading one-row-at-a-time: +# ~~~~~~~~~~~~~~~~~~~~~~~~~~ +# + + +###################################################################### +# Rather than using the full-table methods, if you want to read it +# row-by-row (e.g. to maintain compatibility with an existing event loop), +# you can use a ``TableReader`` instance. +# +# The advantage here is that units and other metadata are retained and +# re-applied +# + + +def read(mode): + + print("reading mode {}".format(mode)) + + with HDF5TableReader("container.h5", mode=mode) as h5_table: + + for group_name in ["data_0/", "data_1/"]: + + group_name = "/{}table".format(group_name) + print(group_name) + + for data in h5_table.read(group_name, VariousTypesContainer): + + print(data.as_dict()) + + +###################################################################### +read("r") + +###################################################################### +read("r+") + +###################################################################### +read("a") + +###################################################################### +read("w") diff --git a/examples/tutorials/README.rst b/examples/tutorials/README.rst new file mode 100644 index 00000000000..c1c31978d8e --- /dev/null +++ b/examples/tutorials/README.rst @@ -0,0 +1,6 @@ +.. _tutorials: + +Tutorials +--------- + +This gallery contains different tutorials of different use cases for ctapipe. diff --git a/examples/tutorials/calibrated_data_exploration.py b/examples/tutorials/calibrated_data_exploration.py new file mode 100644 index 00000000000..40d9d4bb202 --- /dev/null +++ b/examples/tutorials/calibrated_data_exploration.py @@ -0,0 +1,213 @@ +""" +Explore Calibrated Data +======================= + +""" + +import numpy as np +from astropy import units as u +from matplotlib import pyplot as plt + +import ctapipe +from ctapipe.calib import CameraCalibrator +from ctapipe.image import hillas_parameters, tailcuts_clean +from ctapipe.io import EventSource +from ctapipe.utils.datasets import get_dataset_path +from ctapipe.visualization import ArrayDisplay, CameraDisplay + +# %matplotlib inline +plt.style.use("ggplot") + +print(ctapipe.__version__) +print(ctapipe.__file__) + + +###################################################################### +# Let’s first open a raw event file and get an event out of it: +# + +filename = get_dataset_path("gamma_prod5.simtel.zst") +source = EventSource(filename, max_events=2) + +for event in source: + print(event.index.event_id) + +###################################################################### +filename + +###################################################################### +source + +###################################################################### +event + +###################################################################### +print(event.r1) + + +###################################################################### +# Perform basic calibration: +# -------------------------- +# +# Here we will use a ``CameraCalibrator`` which is just a simple wrapper +# that runs the three calibraraton and trace-integration phases of the +# pipeline, taking the data from levels: +# +# **R0** → **R1** → **DL0** → **DL1** +# +# You could of course do these each separately, by using the classes +# ``R1Calibrator``, ``DL0Reducer``, and ``DL1Calibrator``. Note that we +# have not specified any configuration to the ``CameraCalibrator``, so it +# will be using the default algorithms and thresholds, other than +# specifying that the product is a “HESSIOR1Calibrator” (hopefully in the +# near future that will be automatic). +# + + +calib = CameraCalibrator(subarray=source.subarray) +calib(event) + + +###################################################################### +# Now the *r1*, *dl0* and *dl1* containers are filled in the event +# +# - **r1.tel[x]**: contains the “r1-calibrated” waveforms, after +# gain-selection, pedestal subtraciton, and gain-correction +# - **dl0.tel[x]**: is the same but with optional data volume reduction +# (some pixels not filled), in this case this is not performed by +# default, so it is the same as r1 +# - **dl1.tel[x]**: contains the (possibly re-calibrated) waveforms as +# dl0, but also the time-integrated *image* that has been calculated +# using a ``ImageExtractor`` (a ``NeighborPeakWindowSum`` by default) +# + +for tel_id in event.dl1.tel: + print("TEL{:03}: {}".format(tel_id, source.subarray.tel[tel_id])) + print(" - r0 wave shape : {}".format(event.r0.tel[tel_id].waveform.shape)) + print(" - r1 wave shape : {}".format(event.r1.tel[tel_id].waveform.shape)) + print(" - dl1 image shape : {}".format(event.dl1.tel[tel_id].image.shape)) + + +###################################################################### +# Some image processing: +# ---------------------- +# +# Let’s look at the image +# + + +tel_id = sorted(event.r1.tel.keys())[1] +sub = source.subarray +geometry = sub.tel[tel_id].camera.geometry +image = event.dl1.tel[tel_id].image + +###################################################################### +disp = CameraDisplay(geometry, image=image) + +###################################################################### +mask = tailcuts_clean( + geometry, + image, + picture_thresh=10, + boundary_thresh=5, + min_number_picture_neighbors=2, +) +cleaned = image.copy() +cleaned[~mask] = 0 +disp = CameraDisplay(geometry, image=cleaned) + +###################################################################### +params = hillas_parameters(geometry, cleaned) +print(params) +params + +###################################################################### +params = hillas_parameters(geometry, cleaned) + +plt.figure(figsize=(10, 10)) +disp = CameraDisplay(geometry, image=image) +disp.add_colorbar() +disp.overlay_moments(params, color="red", lw=3) +disp.highlight_pixels(mask, color="white", alpha=0.3, linewidth=2) + +plt.xlim(params.x.to_value(u.m) - 0.5, params.x.to_value(u.m) + 0.5) +plt.ylim(params.y.to_value(u.m) - 0.5, params.y.to_value(u.m) + 0.5) + +###################################################################### +source.metadata + + +###################################################################### +# More complex image processing: +# ------------------------------ +# +# Let’s now explore how stereo reconstruction works. +# +# first, look at a summed image from multiple telescopes +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# For this, we want to use a ``CameraDisplay`` again, but since we can’t +# sum and display images with different cameras, we’ll just sub-select +# images from a particular camera type +# +# These are the telescopes that are in this event: +# + +tels_in_event = set( + event.dl1.tel.keys() +) # use a set here, so we can intersect it later +tels_in_event + +###################################################################### +cam_ids = set(sub.get_tel_ids_for_type("MST_MST_NectarCam")) +cam_ids + +###################################################################### +cams_in_event = tels_in_event.intersection(cam_ids) +first_tel_id = list(cams_in_event)[0] +tel = sub.tel[first_tel_id] +print("{}s in event: {}".format(tel, cams_in_event)) + + +###################################################################### +# Now, let’s sum those images: +# + +image_sum = np.zeros_like( + tel.camera.geometry.pix_x.value +) # just make an array of 0's in the same shape as the camera + +for tel_id in cams_in_event: + image_sum += event.dl1.tel[tel_id].image + + +###################################################################### +# And finally display the sum of those images +# + +plt.figure(figsize=(8, 8)) + +disp = CameraDisplay(tel.camera.geometry, image=image_sum) +disp.overlay_moments(params, with_label=False) +plt.title("Sum of {}x {}".format(len(cams_in_event), tel)) + + +###################################################################### +# let’s also show which telescopes those were. Note that currently +# ArrayDisplay’s value field is a vector by ``tel_index``, not ``tel_id``, +# so we have to convert to a tel_index. (this may change in a future +# version to be more user-friendly) +# + + +nectarcam_subarray = sub.select_subarray(cam_ids, name="NectarCam") + +hit_pattern = np.zeros(shape=nectarcam_subarray.n_tels) +hit_pattern[[nectarcam_subarray.tel_indices[x] for x in cams_in_event]] = 100 + +plt.set_cmap(plt.cm.Accent) +plt.figure(figsize=(8, 8)) + +ad = ArrayDisplay(nectarcam_subarray) +ad.values = hit_pattern +ad.add_labels() diff --git a/examples/tutorials/coordinates_example.py b/examples/tutorials/coordinates_example.py new file mode 100644 index 00000000000..3d20aca5735 --- /dev/null +++ b/examples/tutorials/coordinates_example.py @@ -0,0 +1,406 @@ +""" +Coordinates usage in ctapipe +============================ + +""" + +import copy + +import astropy.units as u +import matplotlib.pyplot as plt +from astropy.coordinates import AltAz, EarthLocation, SkyCoord +from astropy.time import Time + +from ctapipe.coordinates import CameraFrame, NominalFrame, TelescopeFrame +from ctapipe.instrument import SubarrayDescription +from ctapipe.io import EventSource +from ctapipe.utils import get_dataset_path +from ctapipe.visualization import CameraDisplay + +# %matplotlib inline + + +# make plots and fonts larger +plt.rcParams["figure.figsize"] = (12, 8) +plt.rcParams["font.size"] = 16 + + +###################################################################### +# Open test dataset +# ----------------- +# + +filename = get_dataset_path("gamma_prod5.simtel.zst") +source = EventSource(filename) + +events = [copy.deepcopy(event) for event in source] +event = events[4] + +layout = set(source.subarray.tel_ids) + + +###################################################################### +# Choose event with LST +# ~~~~~~~~~~~~~~~~~~~~~ +# + + +###################################################################### +# This ensures that the telescope is not “parked” (as it would be in an +# event where it is not triggered) but is actually pointing to a source. +# + +print(f"Telescope with data: {event.r1.tel.keys()}") +tel_id = 3 + + +###################################################################### +# AltAz +# ----- +# +# See `Astropy Docs on +# AltAz `__. +# +# Pointing direction of telescopes or the origin of a simulated shower are +# described in the ``AltAz`` frame. This is a local, angular coordinate +# frame, with angles ``altitude`` and ``azimuth``. Altitude is the +# measured from the Horizon (0°) to the Zenith (90°). For the azimuth, +# there are different conventions. In Astropy und thus ctapipe, Azimuth is +# oriented East of North (i.e., N=0°, E=90°). +# + + +obstime = Time("2013-11-01T03:00") +location = EarthLocation.of_site("Roque de los Muchachos") + +altaz = AltAz(location=location, obstime=obstime) + +array_pointing = SkyCoord( + alt=event.pointing.array_azimuth, + az=event.pointing.array_altitude, + frame=altaz, +) + +print(array_pointing) + + +###################################################################### +# CameraFrame +# ----------- +# +# 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. +# +# **Typical usage**: Position of pixels in the focal plane. +# + +geometry = source.subarray.tel[tel_id].camera.geometry +pix_x = geometry.pix_x +pix_y = geometry.pix_y +focal_length = source.subarray.tel[tel_id].optics.equivalent_focal_length + +###################################################################### +telescope_pointing = SkyCoord( + alt=event.pointing.tel[tel_id].altitude, + az=event.pointing.tel[tel_id].azimuth, + frame=altaz, +) + +camera_frame = CameraFrame( + focal_length=focal_length, + rotation=0 * u.deg, + telescope_pointing=telescope_pointing, +) + +cam_coords = SkyCoord(x=pix_x, y=pix_y, frame=camera_frame) + +print(cam_coords) + +###################################################################### +plt.scatter(cam_coords.x, cam_coords.y) +plt.title(f"Camera type: {geometry.name}") +plt.xlabel(f"x / {cam_coords.x.unit}") +plt.ylabel(f"y / {cam_coords.y.unit}") +plt.axis("square") + + +###################################################################### +# The implementation of the coordinate system with astropy makes it easier +# to use time of the observation and location of the observing site, to +# understand, for example which stars are visible during a certain night +# and how they might be visible in the camera. +# + + +location = EarthLocation.of_site("Roque de los Muchachos") +obstime = Time("2018-11-01T04:00") + +crab = SkyCoord.from_name("crab nebula") + +altaz = AltAz(location=location, obstime=obstime) + +pointing = crab.transform_to(altaz) + +camera_frame = CameraFrame( + telescope_pointing=pointing, + focal_length=focal_length, + obstime=obstime, + location=location, +) + + +subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst") +cam = subarray.tel[1].camera.geometry +fig, ax = plt.subplots() +display = CameraDisplay(cam, ax=ax) + +ax.set_title( + f"La Palma, {obstime}, az={pointing.az.deg:.1f}°, zenith={pointing.zen.deg:.1f}°, camera={geometry.name}" +) + +for i, name in enumerate(["crab nebula", "o tau", "zet tau"]): + star = SkyCoord.from_name(name) + star_cam = star.transform_to(camera_frame) + + x = star_cam.x.to_value(u.m) + y = star_cam.y.to_value(u.m) + + ax.plot(x, y, marker="*", color=f"C{i}") + ax.annotate( + name, + xy=(x, y), + xytext=(5, 5), + textcoords="offset points", + color=f"C{i}", + ) + +plt.show() + + +###################################################################### +# TelescopeFrame +# -------------- +# +# 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 +# ``SkyOffsetFrame``. +# +# The axis of the telescope frame, ``fov_lon`` and ``fov_lat``, are +# aligned with the horizontal system’s azimuth and altitude respectively. +# +# Pointing corrections should applied to the transformation between this +# frame and the camera frame. +# + +telescope_frame = TelescopeFrame( + telescope_pointing=pointing, + obstime=pointing.obstime, + location=pointing.location, +) +telescope_coords = cam_coords.transform_to(telescope_frame) + +###################################################################### +wrap_angle = telescope_pointing.az + 180 * u.deg + +plt.axis("equal") +plt.scatter( + telescope_coords.fov_lon.deg, telescope_coords.fov_lat.deg, alpha=0.2, color="gray" +) + + +for i, name in enumerate(["crab nebula", "o tau", "zet tau"]): + star = SkyCoord.from_name(name) + star_tel = star.transform_to(telescope_frame) + + plt.plot(star_tel.fov_lon.deg, star_tel.fov_lat.deg, "*", ms=10) + plt.annotate( + name, + xy=(star_tel.fov_lon.deg, star_tel.fov_lat.deg), + xytext=(5, 5), + textcoords="offset points", + color=f"C{i}", + ) + +plt.xlabel("fov_lon / {}".format(telescope_coords.altaz.az.unit)) +plt.ylabel("fov_lat / {}".format(telescope_coords.altaz.alt.unit)) + + +###################################################################### +# NominalFrame +# ------------ +# +# 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 ``SkyOffsetFrame`` If the telescopes +# are in divergent pointing, this ``Frame`` can be used to transform to a +# common system. - 2D reconstruction (``HillasIntersector``) is performed +# in this frame - 3D reconstruction (``HillasReconstructor``) doesn’t need +# this frame +# + + +###################################################################### +# Let’s play a bit with 3 LSTs with divergent pointing +# + +location = EarthLocation.of_site("Roque de los Muchachos") +obstime = Time("2018-11-01T02:00") +altaz = AltAz(location=location, obstime=obstime) + +crab = SkyCoord.from_name("crab nebula") + +# let's observe crab +array_pointing = crab.transform_to(altaz) + + +# let the telescopes point to different positions +alt_offsets = u.Quantity([1, -1, -1], u.deg) +az_offsets = u.Quantity([0, -2, +2], u.deg) + + +tel_pointings = SkyCoord( + alt=array_pointing.alt + alt_offsets, + az=array_pointing.az + az_offsets, + frame=altaz, +) + +camera_frames = CameraFrame( + telescope_pointing=tel_pointings, # multiple pointings, so we get multiple frames + focal_length=focal_length, + obstime=obstime, + location=location, +) + +nom_frame = NominalFrame(origin=array_pointing, obstime=obstime, location=location) + +###################################################################### +fig, ax = plt.subplots(figsize=(15, 10)) +ax.set_aspect(1) + +for i in range(3): + cam_coord = SkyCoord(x=pix_x, y=pix_y, frame=camera_frames[i]) + nom_coord = cam_coord.transform_to(nom_frame) + + ax.scatter( + x=nom_coord.fov_lon.deg, + y=nom_coord.fov_lat.deg, + label=f"Telescope {i + 1}", + s=30, + alpha=0.15, + ) + + +for i, name in enumerate(["Crab", "o tau", "zet tau"]): + s = SkyCoord.from_name(name) + s_nom = s.transform_to(nom_frame) + ax.plot( + s_nom.fov_lon.deg, + s_nom.fov_lat.deg, + "*", + ms=10, + ) + ax.annotate( + name, + xy=(s_nom.fov_lon.deg, s_nom.fov_lat.deg), + xytext=(5, 5), + textcoords="offset points", + color=f"C{i}", + ) + + +ax.set_xlabel("fov_lon / deg") +ax.set_ylabel("fov_lat / deg") + +ax.legend() +plt.show() + + +###################################################################### +# GroundFrame +# ----------- +# +# Ground coordinate frame. The ground coordinate frame is a simple +# cartesian frame describing the 3 dimensional position of objects +# compared to the array ground level in relation to the nomial centre of +# the array. Typically this frame will be used for describing the position +# on telescopes and equipment +# +# **Typical usage**: positions of telescopes on the ground (x, y, z) +# + +source.subarray.peek() + + +###################################################################### +# In case a layout is selected, the following line will produce a +# different output from the picture above. +# + +source.subarray.select_subarray(layout, name="Prod3b layout").peek() + + +###################################################################### +# .. figure:: ground_frame.png +# :alt: Ground Frame +# +# Ground Frame + + +###################################################################### +# In this image all the telescope from the ``gamma_test.simtel.gz`` file +# are plotted as spheres in the GroundFrame. +# + + +###################################################################### +# TiltedGroundFrame +# ----------------- +# + + +###################################################################### +# Tilted ground coordinate frame. +# +# The tilted ground coordinate frame is a cartesian system describing the +# 2 dimensional projected positions of objects in a tilted plane described +# by pointing_direction. The plane is rotated along the z_axis by the +# azimuth of the ``pointing_direction`` and then it is inclined with an +# angle equal to the zenith angle of the ``pointing_direction``. +# +# This frame is used for the reconstruction of the shower core position. +# + + +###################################################################### +# .. figure:: tilted_ground_frame.png +# :alt: Tilted Ground Frame +# +# Tilted Ground Frame + + +###################################################################### +# This image picture both the telescopes in the GroundFrame (red) and in +# the TiltedGroundFrame (green) are displayed: in this case since the +# azimuth of the ``pointing_direction`` is 0 degrees, then the plane is +# just tilted according to the zenith angle. +# +# For playing with these and with more 3D models of the telescopes +# themselves, have a look at the +# `CREED_VTK `__ library. +# diff --git a/examples/tutorials/ctapipe_handson.py b/examples/tutorials/ctapipe_handson.py new file mode 100644 index 00000000000..455cf88d9ba --- /dev/null +++ b/examples/tutorials/ctapipe_handson.py @@ -0,0 +1,298 @@ +""" +Getting Started with ctapipe +============================ + +This hands-on was presented at the Paris CTA Consoritum meeting (K. +Kosack) + +""" + + +###################################################################### +# Part 1: load and loop over data +# ------------------------------- +# + +import glob + +import numpy as np +import pandas as pd +from ipywidgets import interact +from matplotlib import pyplot as plt + +from ctapipe import utils +from ctapipe.calib import CameraCalibrator +from ctapipe.image import hillas_parameters, tailcuts_clean +from ctapipe.io import EventSource, HDF5TableWriter +from ctapipe.visualization import CameraDisplay + +# %matplotlib inline + +###################################################################### +path = utils.get_dataset_path("gamma_prod5.simtel.zst") + +###################################################################### +source = EventSource(path, max_events=5) + +for event in source: + print(event.count, event.index.event_id, event.simulation.shower.energy) + +###################################################################### +event + +###################################################################### +event.r1 + +###################################################################### +for event in EventSource(path, max_events=5): + print(event.count, event.r1.tel.keys()) + +###################################################################### +event.r0.tel[3] + +###################################################################### +r0tel = event.r0.tel[3] + +###################################################################### +r0tel.waveform + +###################################################################### +r0tel.waveform.shape + + +###################################################################### +# note that this is (:math:`N_{channels}`, :math:`N_{pixels}`, +# :math:`N_{samples}`) +# + +plt.pcolormesh(r0tel.waveform[0]) + +###################################################################### +brightest_pixel = np.argmax(r0tel.waveform[0].sum(axis=1)) +print(f"pixel {brightest_pixel} has sum {r0tel.waveform[0,1535].sum()}") + +###################################################################### +plt.plot(r0tel.waveform[0, brightest_pixel], label="channel 0 (high-gain)") +plt.plot(r0tel.waveform[1, brightest_pixel], label="channel 1 (low-gain)") +plt.legend() + + +###################################################################### +@interact +def view_waveform(chan=0, pix_id=brightest_pixel): + plt.plot(r0tel.waveform[chan, pix_id]) + + +###################################################################### +# try making this compare 2 waveforms +# + + +###################################################################### +# Part 2: Explore the instrument description +# ------------------------------------------ +# +# This is all well and good, but we don’t really know what camera or +# telescope this is… how do we get instrumental description info? +# +# Currently this is returned *inside* the event (it will soon change to be +# separate in next version or so) +# + +subarray = source.subarray + +###################################################################### +subarray + +###################################################################### +subarray.peek() + +###################################################################### +subarray.to_table() + +###################################################################### +subarray.tel[2] + +###################################################################### +subarray.tel[2].camera + +###################################################################### +subarray.tel[2].optics + +###################################################################### +tel = subarray.tel[2] + +###################################################################### +tel.camera + +###################################################################### +tel.optics + +###################################################################### +tel.camera.geometry.pix_x + +###################################################################### +tel.camera.geometry.to_table() + +###################################################################### +tel.optics.mirror_area + + +###################################################################### +disp = CameraDisplay(tel.camera.geometry) + +###################################################################### +disp = CameraDisplay(tel.camera.geometry) +disp.image = r0tel.waveform[ + 0, :, 10 +] # display channel 0, sample 0 (try others like 10) + + +###################################################################### +# \*\* aside: \*\* show demo using a CameraDisplay in interactive mode in +# ipython rather than notebook +# + + +###################################################################### +# Part 3: Apply some calibration and trace integration +# ---------------------------------------------------- +# + + +calib = CameraCalibrator(subarray=subarray) + +###################################################################### +for event in EventSource(path, max_events=5): + calib(event) # fills in r1, dl0, and dl1 + print(event.dl1.tel.keys()) + +###################################################################### +event.dl1.tel[3] + +###################################################################### +dl1tel = event.dl1.tel[3] + +###################################################################### +dl1tel.image.shape # note this will be gain-selected in next version, so will be just 1D array of 1855 + +###################################################################### +dl1tel.peak_time + +###################################################################### +CameraDisplay(tel.camera.geometry, image=dl1tel.image) + +###################################################################### +CameraDisplay(tel.camera.geometry, image=dl1tel.peak_time) + + +###################################################################### +# Now for Hillas Parameters +# + + +image = dl1tel.image +mask = tailcuts_clean(tel.camera.geometry, image, picture_thresh=10, boundary_thresh=5) +mask + +###################################################################### +CameraDisplay(tel.camera.geometry, image=mask) + +###################################################################### +cleaned = image.copy() +cleaned[~mask] = 0 + +###################################################################### +disp = CameraDisplay(tel.camera.geometry, image=cleaned) +disp.cmap = plt.cm.coolwarm +disp.add_colorbar() +plt.xlim(0.5, 1.0) +plt.ylim(-1.0, 0.0) + +###################################################################### +params = hillas_parameters(tel.camera.geometry, cleaned) +print(params) + +###################################################################### +disp = CameraDisplay(tel.camera.geometry, image=cleaned) +disp.cmap = plt.cm.coolwarm +disp.add_colorbar() +plt.xlim(0.5, 1.0) +plt.ylim(-1.0, 0.0) +disp.overlay_moments(params, color="white", lw=2) + + +###################################################################### +# Part 4: Let’s put it all together: +# ---------------------------------- +# +# - loop over events, selecting only telescopes of the same type +# (e.g. LST:LSTCam) +# - for each event, apply calibration/trace integration +# - calculate Hillas parameters +# - write out all hillas paremeters to a file that can be loaded with +# Pandas +# + + +###################################################################### +# first let’s select only those telescopes with LST:LSTCam +# + +subarray.telescope_types + +###################################################################### +subarray.get_tel_ids_for_type("LST_LST_LSTCam") + + +###################################################################### +# Now let’s write out program +# + +data = utils.get_dataset_path("gamma_prod5.simtel.zst") +source = EventSource(data) # remove the max_events limit to get more stats + +###################################################################### +for event in source: + calib(event) + + for tel_id, tel_data in event.dl1.tel.items(): + tel = source.subarray.tel[tel_id] + mask = tailcuts_clean(tel.camera.geometry, tel_data.image) + if np.count_nonzero(mask) > 0: + params = hillas_parameters(tel.camera.geometry[mask], tel_data.image[mask]) + + +###################################################################### +with HDF5TableWriter(filename="hillas.h5", group_name="dl1", overwrite=True) as writer: + + source = EventSource(data, allowed_tels=[1, 2, 3, 4], max_events=10) + for event in source: + calib(event) + + for tel_id, tel_data in event.dl1.tel.items(): + tel = source.subarray.tel[tel_id] + mask = tailcuts_clean(tel.camera.geometry, tel_data.image) + params = hillas_parameters(tel.camera.geometry[mask], tel_data.image[mask]) + writer.write("hillas", params) + + +###################################################################### +# We can now load in the file we created and plot it +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +glob.glob("*.h5") + + +###################################################################### +hillas = pd.read_hdf("hillas.h5", key="/dl1/hillas") +hillas + +###################################################################### +_ = hillas.hist(figsize=(8, 8)) + + +###################################################################### +# If you do this yourself, chose a larger file to loop over more events to +# get better statistics +# diff --git a/examples/tutorials/ctapipe_overview.py b/examples/tutorials/ctapipe_overview.py new file mode 100644 index 00000000000..e96cea70f45 --- /dev/null +++ b/examples/tutorials/ctapipe_overview.py @@ -0,0 +1,653 @@ +""" +Analyzing Events Using ctapipe +============================== + +""" + + +###################################################################### +# Initially presented @ LST Analysis Bootcamp +# in Padova, 26.11.2018 +# by Maximilian Nöthe (@maxnoe) & Kai A. Brügge (@mackaiver) + + +import tempfile +import timeit +from copy import deepcopy + +import astropy.units as u +import matplotlib.pyplot as plt +import numpy as np +from astropy.coordinates import AltAz +from astropy.coordinates.angle_utilities import angular_separation +from matplotlib.colors import ListedColormap +from scipy.sparse.csgraph import connected_components +from traitlets.config import Config + +from ctapipe.calib import CameraCalibrator +from ctapipe.image import ( + ImageProcessor, + camera_to_shower_coordinates, + concentration_parameters, + hillas_parameters, + leakage_parameters, + number_of_islands, + timing_parameters, + toymodel, +) +from ctapipe.image.cleaning import tailcuts_clean +from ctapipe.io import DataWriter, EventSource, TableLoader +from ctapipe.reco import ShowerProcessor +from ctapipe.utils.datasets import get_dataset_path +from ctapipe.visualization import ArrayDisplay, CameraDisplay + +# %matplotlib inline + +###################################################################### +plt.rcParams["figure.figsize"] = (12, 8) +plt.rcParams["font.size"] = 14 +plt.rcParams["figure.figsize"] + + +###################################################################### +# General Information +# ------------------- +# + + +###################################################################### +# Design +# ~~~~~~ +# +# - DL0 → DL3 analysis +# +# - Currently some R0 → DL2 code to be able to analyze simtel files +# +# - ctapipe is built upon the Scientific Python Stack, core dependencies +# are +# +# - numpy +# - scipy +# - astropy +# - numba +# + + +###################################################################### +# Developement +# ~~~~~~~~~~~~ +# +# - ctapipe is developed as Open Source Software (BSD 3-Clause License) +# at https://github.com/cta-observatory/ctapipe +# +# - We use the “Github-Workflow”: +# +# - Few people (e.g. @kosack, @maxnoe) have write access to the main +# repository +# - Contributors fork the main repository and work on branches +# - Pull Requests are merged after Code Review and automatic execution +# of the test suite +# +# - Early developement stage ⇒ backwards-incompatible API changes might +# and will happen +# + + +###################################################################### +# What’s there? +# ~~~~~~~~~~~~~ +# +# - Reading simtel simulation files +# - Simple calibration, cleaning and feature extraction functions +# - Camera and Array plotting +# - Coordinate frames and transformations +# - Stereo-reconstruction using line intersections +# + + +###################################################################### +# What’s still missing? +# ~~~~~~~~~~~~~~~~~~~~~ +# +# - Good integration with machine learning techniques +# - IRF calculation +# - Documentation, e.g. formal definitions of coordinate frames +# + + +###################################################################### +# What can you do? +# ~~~~~~~~~~~~~~~~ +# +# - Report issues +# +# - Hard to get started? Tell us where you are stuck +# - Tell user stories +# - Missing features +# +# - Start contributing +# +# - ctapipe needs more workpower +# - Implement new reconstruction features +# + + +###################################################################### +# A simple hillas analysis +# ------------------------ +# + + +###################################################################### +# Reading in simtel files +# ~~~~~~~~~~~~~~~~~~~~~~~ +# + + +input_url = get_dataset_path("gamma_prod5.simtel.zst") + +# EventSource() automatically detects what kind of file we are giving it, +# if already supported by ctapipe +source = EventSource(input_url, max_events=5) + +print(type(source)) + +###################################################################### +for event in source: + print( + "Id: {}, E = {:1.3f}, Telescopes: {}".format( + event.count, event.simulation.shower.energy, len(event.r0.tel) + ) + ) + + +###################################################################### +# Each event is a ``DataContainer`` holding several ``Field``\ s of data, +# which can be containers or just numbers. Let’s look a one event: +# + +event + +###################################################################### +source.subarray.camera_types + +###################################################################### +len(event.r0.tel), len(event.r1.tel) + + +###################################################################### +# Data calibration +# ~~~~~~~~~~~~~~~~ +# +# The ``CameraCalibrator`` calibrates the event (obtaining the ``dl1`` +# images). +# + + +calibrator = CameraCalibrator(subarray=source.subarray) + +###################################################################### +calibrator(event) + + +###################################################################### +# Event displays +# ~~~~~~~~~~~~~~ +# +# Let’s use ctapipe’s plotting facilities to plot the telescope images +# + +event.dl1.tel.keys() + +###################################################################### +tel_id = 130 + +###################################################################### +geometry = source.subarray.tel[tel_id].camera.geometry +dl1 = event.dl1.tel[tel_id] + +geometry, dl1 + +###################################################################### +dl1.image + + +###################################################################### +display = CameraDisplay(geometry) + +# right now, there might be one image per gain channel. +# This will change as soon as +display.image = dl1.image +display.add_colorbar() + + +###################################################################### +# Image Cleaning +# ~~~~~~~~~~~~~~ +# + + +# unoptimized cleaning levels +cleaning_level = { + "CHEC": (2, 4, 2), + "LSTCam": (3.5, 7, 2), + "FlashCam": (3.5, 7, 2), + "NectarCam": (4, 8, 2), +} + +###################################################################### +boundary, picture, min_neighbors = cleaning_level[geometry.name] + +clean = tailcuts_clean( + geometry, + dl1.image, + boundary_thresh=boundary, + picture_thresh=picture, + min_number_picture_neighbors=min_neighbors, +) + +###################################################################### +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5)) + +d1 = CameraDisplay(geometry, ax=ax1) +d2 = CameraDisplay(geometry, ax=ax2) + +ax1.set_title("Image") +d1.image = dl1.image +d1.add_colorbar(ax=ax1) + +ax2.set_title("Pulse Time") +d2.image = dl1.peak_time - np.average(dl1.peak_time, weights=dl1.image) +d2.cmap = "RdBu_r" +d2.add_colorbar(ax=ax2) +d2.set_limits_minmax(-20, 20) + +d1.highlight_pixels(clean, color="red", linewidth=1) + + +###################################################################### +# Image Parameters +# ~~~~~~~~~~~~~~~~ +# + + +hillas = hillas_parameters(geometry[clean], dl1.image[clean]) + +print(hillas) + +###################################################################### +display = CameraDisplay(geometry) + +# set "unclean" pixels to 0 +cleaned = dl1.image.copy() +cleaned[~clean] = 0.0 + +display.image = cleaned +display.add_colorbar() + +display.overlay_moments(hillas, color="xkcd:red") + +###################################################################### +timing = timing_parameters(geometry, dl1.image, dl1.peak_time, hillas, clean) + +print(timing) + +###################################################################### +long, trans = camera_to_shower_coordinates( + geometry.pix_x, geometry.pix_y, hillas.x, hillas.y, hillas.psi +) + +plt.plot(long[clean], dl1.peak_time[clean], "o") +plt.plot(long[clean], timing.slope * long[clean] + timing.intercept) + +###################################################################### +leakage = leakage_parameters(geometry, dl1.image, clean) +print(leakage) + +###################################################################### +disp = CameraDisplay(geometry) +disp.image = dl1.image +disp.highlight_pixels(geometry.get_border_pixel_mask(1), linewidth=2, color="xkcd:red") + +###################################################################### +n_islands, island_id = number_of_islands(geometry, clean) + +print(n_islands) + +###################################################################### +conc = concentration_parameters(geometry, dl1.image, hillas) +print(conc) + + +###################################################################### +# Putting it all together / Stereo reconstruction +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# All these steps are now unified in several components configurable +# through the config system, mainly: +# +# - CameraCalibrator for DL0 → DL1 (Images) +# - ImageProcessor for DL1 (Images) → DL1 (Parameters) +# - ShowerProcessor for stereo reconstruction of the shower geometry +# - DataWriter for writing data into HDF5 +# +# A command line tool doing these steps and writing out data in HDF5 +# format is available as ``ctapipe-process`` +# + + +image_processor_config = Config( + { + "ImageProcessor": { + "image_cleaner_type": "TailcutsImageCleaner", + "TailcutsImageCleaner": { + "picture_threshold_pe": [ + ("type", "LST_LST_LSTCam", 7.5), + ("type", "MST_MST_FlashCam", 8), + ("type", "MST_MST_NectarCam", 8), + ("type", "SST_ASTRI_CHEC", 7), + ], + "boundary_threshold_pe": [ + ("type", "LST_LST_LSTCam", 5), + ("type", "MST_MST_FlashCam", 4), + ("type", "MST_MST_NectarCam", 4), + ("type", "SST_ASTRI_CHEC", 4), + ], + }, + } + } +) + +input_url = get_dataset_path("gamma_prod5.simtel.zst") +source = EventSource(input_url) + +calibrator = CameraCalibrator(subarray=source.subarray) +image_processor = ImageProcessor( + subarray=source.subarray, config=image_processor_config +) +shower_processor = ShowerProcessor(subarray=source.subarray) +horizon_frame = AltAz() + +f = tempfile.NamedTemporaryFile(suffix=".hdf5") + +with DataWriter( + source, output_path=f.name, overwrite=True, write_showers=True +) as writer: + + for event in source: + energy = event.simulation.shower.energy + n_telescopes_r1 = len(event.r1.tel) + event_id = event.index.event_id + print(f"Id: {event_id}, E = {energy:1.3f}, Telescopes (R1): {n_telescopes_r1}") + + calibrator(event) + image_processor(event) + shower_processor(event) + + stereo = event.dl2.stereo.geometry["HillasReconstructor"] + if stereo.is_valid: + print(" Alt: {:.2f}°".format(stereo.alt.deg)) + print(" Az: {:.2f}°".format(stereo.az.deg)) + print(" Hmax: {:.0f}".format(stereo.h_max)) + print(" CoreX: {:.1f}".format(stereo.core_x)) + print(" CoreY: {:.1f}".format(stereo.core_y)) + print(" Multiplicity: {:d}".format(len(stereo.telescopes))) + + # save a nice event for plotting later + if event.count == 3: + plotting_event = deepcopy(event) + + writer(event) + + +###################################################################### +loader = TableLoader(f.name, load_dl2=True, load_simulated=True) + +events = loader.read_subarray_events() + +###################################################################### +theta = angular_separation( + events["HillasReconstructor_az"].quantity, + events["HillasReconstructor_alt"].quantity, + events["true_az"].quantity, + events["true_alt"].quantity, +) + +plt.hist(theta.to_value(u.deg) ** 2, bins=25, range=[0, 0.3]) +plt.xlabel(r"$\theta² / deg²$") +None + + +###################################################################### +# ArrayDisplay +# ------------ +# + + +angle_offset = plotting_event.pointing.array_azimuth + +plotting_hillas = { + tel_id: dl1.parameters.hillas for tel_id, dl1 in plotting_event.dl1.tel.items() +} + +plotting_core = { + tel_id: dl1.parameters.core.psi for tel_id, dl1 in plotting_event.dl1.tel.items() +} + + +disp = ArrayDisplay(source.subarray) + +disp.set_line_hillas(plotting_hillas, plotting_core, 500) + +plt.scatter( + plotting_event.simulation.shower.core_x, + plotting_event.simulation.shower.core_y, + s=200, + c="k", + marker="x", + label="True Impact", +) +plt.scatter( + plotting_event.dl2.stereo.geometry["HillasReconstructor"].core_x, + plotting_event.dl2.stereo.geometry["HillasReconstructor"].core_y, + s=200, + c="r", + marker="x", + label="Estimated Impact", +) + +plt.legend() +# plt.xlim(-400, 400) +# plt.ylim(-400, 400) + + +###################################################################### +# Reading the LST dl1 data +# ~~~~~~~~~~~~~~~~~~~~~~~~ +# + +loader = TableLoader(f.name, load_simulated=True, load_dl1_parameters=True) + +dl1_table = loader.read_telescope_events(["LST_LST_LSTCam"]) + +###################################################################### +plt.scatter( + np.log10(dl1_table["true_energy"].quantity / u.TeV), + np.log10(dl1_table["hillas_intensity"]), +) +plt.xlabel("log10(E / TeV)") +plt.ylabel("log10(intensity)") +None + + +###################################################################### +# Isn’t python slow? +# ------------------ +# +# - Many of you might have heard: “Python is slow”. +# - That’s trueish. +# - All python objects are classes living on the heap, even integers. +# - Looping over lots of “primitives” is quite slow compared to other +# languages. +# +# | ⇒ Vectorize as much as possible using numpy +# | ⇒ Use existing interfaces to fast C / C++ / Fortran code +# | ⇒ Optimize using numba +# +# **But: “Premature Optimization is the root of all evil” — Donald Knuth** +# +# So profile to find exactly what is slow. +# +# Why use python then? +# ~~~~~~~~~~~~~~~~~~~~ +# +# - Python works very well as *glue* for libraries of all kinds of +# languages +# - Python has a rich ecosystem for data science, physics, algorithms, +# astronomy +# +# Example: Number of Islands +# ~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Find all groups of pixels, that survived the cleaning +# + + +geometry = loader.subarray.tel[1].camera.geometry + + +###################################################################### +# Let’s create a toy images with several islands; +# + +np.random.seed(42) + +image = np.zeros(geometry.n_pixels) + + +for i in range(9): + + model = toymodel.Gaussian( + x=np.random.uniform(-0.8, 0.8) * u.m, + y=np.random.uniform(-0.8, 0.8) * u.m, + width=np.random.uniform(0.05, 0.075) * u.m, + length=np.random.uniform(0.1, 0.15) * u.m, + psi=np.random.uniform(0, 2 * np.pi) * u.rad, + ) + + new_image, sig, bg = model.generate_image( + geometry, intensity=np.random.uniform(1000, 3000), nsb_level_pe=5 + ) + image += new_image + +###################################################################### +clean = tailcuts_clean( + geometry, + image, + picture_thresh=10, + boundary_thresh=5, + min_number_picture_neighbors=2, +) + +###################################################################### +disp = CameraDisplay(geometry) +disp.image = image +disp.highlight_pixels(clean, color="xkcd:red", linewidth=1.5) +disp.add_colorbar() + + +###################################################################### +def num_islands_python(camera, clean): + """A breadth first search to find connected islands of neighboring pixels in the cleaning set""" + + # the camera geometry has a [n_pixel, n_pixel] boolean array + # that is True where two pixels are neighbors + neighbors = camera.neighbor_matrix + + island_ids = np.zeros(camera.n_pixels) + current_island = 0 + + # a set to remember which pixels we already visited + visited = set() + + # go only through the pixels, that survived cleaning + for pix_id in np.where(clean)[0]: + if pix_id not in visited: + # remember that we already checked this pixel + visited.add(pix_id) + + # if we land in the outer loop again, we found a new island + current_island += 1 + island_ids[pix_id] = current_island + + # now check all neighbors of the current pixel recursively + to_check = set(np.where(neighbors[pix_id] & clean)[0]) + while to_check: + pix_id = to_check.pop() + + if pix_id not in visited: + visited.add(pix_id) + island_ids[pix_id] = current_island + + to_check.update(np.where(neighbors[pix_id] & clean)[0]) + + n_islands = current_island + return n_islands, island_ids + + +###################################################################### +n_islands, island_ids = num_islands_python(geometry, clean) + + +###################################################################### +cmap = plt.get_cmap("Paired") +cmap = ListedColormap(cmap.colors[:n_islands]) +cmap.set_under("k") + +disp = CameraDisplay(geometry) +disp.image = island_ids +disp.cmap = cmap +disp.set_limits_minmax(0.5, n_islands + 0.5) +disp.add_colorbar() + +###################################################################### +timeit.timeit(lambda: num_islands_python(geometry, clean), number=1000) / 1000 + + +###################################################################### +def num_islands_scipy(geometry, clean): + neighbors = geometry.neighbor_matrix_sparse + + clean_neighbors = neighbors[clean][:, clean] + num_islands, labels = connected_components(clean_neighbors, directed=False) + + island_ids = np.zeros(geometry.n_pixels) + island_ids[clean] = labels + 1 + + return num_islands, island_ids + + +###################################################################### +n_islands_s, island_ids_s = num_islands_scipy(geometry, clean) + +###################################################################### +disp = CameraDisplay(geometry) +disp.image = island_ids_s +disp.cmap = cmap +disp.set_limits_minmax(0.5, n_islands_s + 0.5) +disp.add_colorbar() + +###################################################################### +timeit.timeit(lambda: num_islands_scipy(geometry, clean), number=10000) / 10000 + +###################################################################### +# **A lot less code, and a factor 3 speed improvement** +# + + +###################################################################### +# Finally, current ctapipe implementation is using numba: +# + +###################################################################### +timeit.timeit(lambda: number_of_islands(geometry, clean), number=100000) / 100000 diff --git a/docs/user-guide/tutorials/ground_frame.png b/examples/tutorials/ground_frame.png similarity index 100% rename from docs/user-guide/tutorials/ground_frame.png rename to examples/tutorials/ground_frame.png diff --git a/examples/tutorials/raw_data_exploration.py b/examples/tutorials/raw_data_exploration.py new file mode 100644 index 00000000000..7fa6dc3c6b4 --- /dev/null +++ b/examples/tutorials/raw_data_exploration.py @@ -0,0 +1,308 @@ +""" +Exploring Raw Data +================== + +Here are just some very simple examples of going through and inspecting +the raw data, and making some plots using ``ctapipe``. The data explored +here are *raw Monte Carlo* data, which is Data Level “R0” in CTA +terminology (e.g. it is before any processing that would happen inside a +Camera or off-line) + +""" + + +###################################################################### +# Setup: +# + +import numpy as np +from matplotlib import pyplot as plt +from scipy import signal + +from ctapipe.io import EventSource +from ctapipe.utils import get_dataset_path +from ctapipe.visualization import CameraDisplay + +# %matplotlib inline + + +###################################################################### +# To read SimTelArray format data, ctapipe uses the ``pyeventio`` library +# (which is installed automatically along with ctapipe). The following +# lines however will load any data known to ctapipe (multiple +# ``EventSources`` are implemented, and chosen automatically based on the +# type of the input file. +# +# All data access first starts with an ``EventSource``, and here we use a +# helper function ``event_source`` that constructs one. The resulting +# ``source`` object can be iterated over like a list of events. We also +# here use an ``EventSeeker`` which provides random-access to the source +# (by seeking to the given event ID or number) +# + +source = EventSource(get_dataset_path("gamma_prod5.simtel.zst"), max_events=5) + + +###################################################################### +# Explore the contents of an event +# -------------------------------- +# +# note that the R0 level is the raw data that comes out of a camera, and +# also the lowest level of monte-carlo data. +# + +# so we can advance through events one-by-one +event_iterator = iter(source) + +event = next(event_iterator) + + +###################################################################### +# the event is just a class with a bunch of data items in it. You can see +# a more compact represntation via: +# + +event.r0 + + +###################################################################### +# printing the event structure, will currently print the value all items +# under it (so you get a lot of output if you print a high-level +# container): +# + +print(event.simulation.shower) + +###################################################################### +print(event.r0.tel.keys()) + + +###################################################################### +# note that the event has 3 telescopes in it: Let’s try the next one: +# + +event = next(event_iterator) +print(event.r0.tel.keys()) + + +###################################################################### +# now, we have a larger event with many telescopes… Let’s look at one of +# them: +# + +teldata = event.r0.tel[26] +print(teldata) +teldata + + +###################################################################### +# Note that some values are unit quantities (``astropy.units.Quantity``) +# or angular quantities (``astropy.coordinates.Angle``), and you can +# easily maniuplate them: +# + +event.simulation.shower.energy + +###################################################################### +event.simulation.shower.energy.to("GeV") + +###################################################################### +event.simulation.shower.energy.to("J") + +###################################################################### +event.simulation.shower.alt + +###################################################################### +print("Altitude in degrees:", event.simulation.shower.alt.deg) + + +###################################################################### +# Look for signal pixels in a camera +# ---------------------------------- +# +# again, ``event.r0.tel[x]`` contains a data structure for the telescope +# data, with some fields like ``waveform``. +# +# Let’s make a 2D plot of the sample data (sample vs pixel), so we can see +# if we see which pixels contain Cherenkov light signals: +# + +plt.pcolormesh(teldata.waveform[0]) # note the [0] is for channel 0 +plt.colorbar() +plt.xlabel("sample number") +plt.ylabel("Pixel_id") + + +###################################################################### +# Let’s zoom in to see if we can identify the pixels that have the +# Cherenkov signal in them +# + +plt.pcolormesh(teldata.waveform[0]) +plt.colorbar() +plt.ylim(700, 750) +plt.xlabel("sample number") +plt.ylabel("pixel_id") +print("waveform[0] is an array of shape (N_pix,N_slice) =", teldata.waveform[0].shape) + + +###################################################################### +# Now we can really see that some pixels have a signal in them! +# +# Lets look at a 1D plot of pixel 270 in channel 0 and see the signal: +# + +trace = teldata.waveform[0][719] +plt.plot(trace, drawstyle="steps") + + +###################################################################### +# Great! It looks like a *standard Cherenkov signal*! +# +# Let’s take a look at several traces to see if the peaks area aligned: +# + +for pix_id in range(718, 723): + plt.plot( + teldata.waveform[0][pix_id], label="pix {}".format(pix_id), drawstyle="steps" + ) +plt.legend() + + +###################################################################### +# Look at the time trace from a Camera Pixel +# ------------------------------------------ +# +# ``ctapipe.calib.camera`` includes classes for doing automatic trace +# integration with many methods, but before using that, let’s just try to +# do something simple! +# +# Let’s define the integration windows first: By eye, they seem to be +# reaonsable from sample 8 to 13 for signal, and 20 to 29 for pedestal +# (which we define as the sum of all noise: NSB + electronic) +# + +for pix_id in range(718, 723): + plt.plot(teldata.waveform[0][pix_id], "+-") +plt.fill_betweenx([0, 1600], 19, 24, color="red", alpha=0.3, label="Ped window") +plt.fill_betweenx([0, 1600], 5, 9, color="green", alpha=0.3, label="Signal window") +plt.legend() + + +###################################################################### +# Do a very simplisitic trace analysis +# ------------------------------------ +# +# Now, let’s for example calculate a signal and background in a the fixed +# windows we defined for this single event. Note we are ignoring the fact +# that cameras have 2 gains, and just using a single gain (channel 0, +# which is the high-gain channel): +# + +data = teldata.waveform[0] +peds = data[:, 19:24].mean(axis=1) # mean of samples 20 to 29 for all pixels +sums = data[:, 5:9].sum(axis=1) / (13 - 8) # simple sum integration + +###################################################################### +phist = plt.hist(peds, bins=50, range=[0, 150]) +plt.title("Pedestal Distribution of all pixels for a single event") + + +###################################################################### +# let’s now take a look at the pedestal-subtracted sums and a +# pedestal-subtracted signal: +# + +plt.plot(sums - peds) +plt.xlabel("pixel id") +plt.ylabel("Pedestal-subtracted Signal") + + +###################################################################### +# Now, we can clearly see that the signal is centered at 0 where there is +# no Cherenkov light, and we can also clearly see the shower around pixel +# 250. +# + +# we can also subtract the pedestals from the traces themselves, which would be needed to compare peaks properly +for ii in range(270, 280): + plt.plot(data[ii] - peds[ii], drawstyle="steps", label="pix{}".format(ii)) +plt.legend() + + +###################################################################### +# Camera Displays +# --------------- +# +# It’s of course much easier to see the signal if we plot it in 2D with +# correct pixel positions! +# +# note: the instrument data model is not fully implemented, so there is +# not a good way to load all the camera information (right now it is +# hacked into the ``inst`` sub-container that is read from the +# Monte-Carlo file) +# + +camgeom = source.subarray.tel[24].camera.geometry + +###################################################################### +title = "CT24, run {} event {} ped-sub".format(event.index.obs_id, event.index.event_id) +disp = CameraDisplay(camgeom, title=title) +disp.image = sums - peds +disp.cmap = plt.cm.RdBu_r +disp.add_colorbar() +disp.set_limits_percent(95) # autoscale + + +###################################################################### +# It looks like a nice signal! We have plotted our pedestal-subtracted +# trace integral, and see the shower clearly! +# +# Let’s look at all telescopes: +# +# note we plot here the raw signal, since we have not calculated the +# pedestals for each) +# + +for tel in event.r0.tel.keys(): + plt.figure() + camgeom = source.subarray.tel[tel].camera.geometry + title = "CT{}, run {} event {}".format( + tel, event.index.obs_id, event.index.event_id + ) + disp = CameraDisplay(camgeom, title=title) + disp.image = event.r0.tel[tel].waveform[0].sum(axis=1) + disp.cmap = plt.cm.RdBu_r + disp.add_colorbar() + disp.set_limits_percent(95) + + +###################################################################### +# some signal processing… +# ----------------------- +# +# Let’s try to detect the peak using the scipy.signal package: +# https://docs.scipy.org/doc/scipy/reference/signal.html +# + + +pix_ids = np.arange(len(data)) +has_signal = sums > 300 + +widths = np.array( + [ + 8, + ] +) # peak widths to search for (let's fix it at 8 samples, about the width of the peak) +peaks = [signal.find_peaks_cwt(trace, widths) for trace in data[has_signal]] + +for p, s in zip(pix_ids[has_signal], peaks): + print("pix{} has peaks at sample {}".format(p, s)) + plt.plot(data[p], drawstyle="steps-mid") + plt.scatter(np.array(s), data[p, s]) + + +###################################################################### +# clearly the signal needs to be filtered first, or an appropriate wavelet +# used, but the idea is nice +# diff --git a/examples/tutorials/theta_square.py b/examples/tutorials/theta_square.py new file mode 100644 index 00000000000..b7f9694383e --- /dev/null +++ b/examples/tutorials/theta_square.py @@ -0,0 +1,76 @@ +r""" +Make a theta-square plot +======================== + +This is a basic example to analyze some events and make a +:math:`\Theta^2` plot + +""" + +# %matplotlib inline + +import matplotlib.pyplot as plt +import numpy as np +from astropy import units as u +from astropy.coordinates.angle_utilities import angular_separation +from tqdm.auto import tqdm + +from ctapipe.calib import CameraCalibrator +from ctapipe.image import ImageProcessor +from ctapipe.io import EventSource +from ctapipe.reco import ShowerProcessor + +###################################################################### +# Get source events in MC dataset. +# + +source = EventSource( + "dataset://gamma_prod5.simtel.zst", + # allowed_tels={1, 2, 3, 4}, +) + +subarray = source.subarray + +calib = CameraCalibrator(subarray=subarray) +image_processor = ImageProcessor(subarray=subarray) +shower_processor = ShowerProcessor(subarray=subarray) + +off_angles = [] + +for event in tqdm(source): + + # calibrating the event + calib(event) + image_processor(event) + shower_processor(event) + + reco_result = event.dl2.stereo.geometry["HillasReconstructor"] + + # get angular offset between reconstructed shower direction and MC + # generated shower direction + true_shower = event.simulation.shower + off_angle = angular_separation( + true_shower.az, true_shower.alt, reco_result.az, reco_result.alt + ) + + # Appending all estimated off angles + off_angles.append(off_angle.to(u.deg).value) + + +###################################################################### +# calculate theta square for angles which are not nan +# + +off_angles = np.array(off_angles) +thetasquare = off_angles[np.isfinite(off_angles)] ** 2 + + +###################################################################### +# Plot the results +# ---------------- +# + +plt.hist(thetasquare, bins=10, range=[0, 0.4]) +plt.xlabel(r"$\theta^2$ (deg)") +plt.ylabel("# of events") +plt.show() diff --git a/docs/user-guide/tutorials/tilted_ground_frame.png b/examples/tutorials/tilted_ground_frame.png similarity index 100% rename from docs/user-guide/tutorials/tilted_ground_frame.png rename to examples/tutorials/tilted_ground_frame.png diff --git a/examples/visualization/README.rst b/examples/visualization/README.rst new file mode 100644 index 00000000000..ca55d616351 --- /dev/null +++ b/examples/visualization/README.rst @@ -0,0 +1,4 @@ +.. _visualization-examples-gallery: + +Visualization +------------- diff --git a/examples/visualization/array_display.py b/examples/visualization/array_display.py new file mode 100644 index 00000000000..3dff647f7ae --- /dev/null +++ b/examples/visualization/array_display.py @@ -0,0 +1,212 @@ +""" +Array Displays +============== + +Like ``CameraDisplays``, ctapipe provides a way to display information +related to the array on the ground: ``ArrayDisplay`` + +""" + +import matplotlib.pyplot as plt +import numpy as np +from astropy import units as u +from astropy.coordinates import SkyCoord + +from ctapipe.calib import CameraCalibrator +from ctapipe.coordinates import EastingNorthingFrame +from ctapipe.image import ImageProcessor +from ctapipe.instrument import SubarrayDescription +from ctapipe.io import EventSource +from ctapipe.reco import ShowerProcessor +from ctapipe.visualization import ArrayDisplay + +plt.rcParams["figure.figsize"] = (8, 6) + +###################################################################### +tel_ids = list(range(1, 5)) + list(range(5, 20)) # just LSTs + one set of MSTs + +subarray = SubarrayDescription.read( + "dataset://gamma_20deg_0deg_run1___cta-prod5-lapalma_desert-2158m-LaPalma-dark_100evts.simtel.zst" +).select_subarray(tel_ids) + + +###################################################################### +# An array display is created for example in ``subarray.peek()``: +# + +subarray.peek() + + +###################################################################### +# However, you can make one manually with a bit more flexibility: +# + + +###################################################################### +# Constructing an ArrayDisplay +# ---------------------------- +# + +disp = ArrayDisplay(subarray) + + +###################################################################### +# You can specify the Frame you want as long as it is compatible with +# ``GroundFrame``. ``EastingNorthingFrame`` is probably the most useful. +# You can also add telescope labels +# + +disp = ArrayDisplay(subarray, frame=EastingNorthingFrame()) +disp.add_labels() + + +###################################################################### +# Using color to show information +# ------------------------------- +# +# By default the color of the telescope circles correlates to telescope +# type. However, you can use color to convey other information by setting +# the ``values`` attribute, like a trigger pattern +# + +plt.set_cmap("rainbow") # the array display will use the current colormap for values + +ad = ArrayDisplay(subarray) +ad.telescopes.set_linewidth(0) # to turn off the telescope borders + +trigger_pattern = np.zeros(subarray.n_tels) +trigger_pattern[ + [ + 1, + 4, + 5, + 6, + ] +] = 1 +ad.values = trigger_pattern # display certain telescopes in a color +ad.add_labels() + + +###################################################################### +# or for example, you could use color to represent the telescope distance +# to the impact point +# + +shower_impact = SkyCoord(200 * u.m, -200 * u.m, 0 * u.m, frame=EastingNorthingFrame()) + +plt.set_cmap("rainbow") # the array display will use the current colormap for values +ad = ArrayDisplay(subarray) +ad.telescopes.set_linewidth(0) # to turn off the telescope borders +plt.scatter(shower_impact.easting, shower_impact.northing, marker="+", s=200) + +distances = np.hypot( + subarray.tel_coords.cartesian.x - shower_impact.cartesian.x, + subarray.tel_coords.cartesian.y - shower_impact.cartesian.y, +) +ad.values = distances +plt.colorbar(ad.telescopes, label="Distance (m)") + + +###################################################################### +# Overlaying vectors +# ------------------ +# +# For plotting reconstruction quantities, it’s useful to overlay vectors +# on the telescope positions. ``ArrayDisplay`` provides functions: \* +# ``set_vector_uv`` to set by cartesian coordinates from the center of +# each telescope \* ``set_vector_rho_phi`` to set by polar coorinates from +# the center of each telescope \* ``set_vector_hillas`` to set vectors +# from a ``dict[int,HillasParameters]`` mapping tel_id (not index!) to a +# set of parameters. +# + +np.random.seed(0) +phis = np.random.uniform(0, 180.0, size=subarray.n_tels) * u.deg +rhos = np.ones(subarray.n_tels) * 200 * u.m + + +ad = ArrayDisplay(subarray, frame=EastingNorthingFrame(), tel_scale=2) +ad.set_vector_rho_phi(rho=rhos, phi=phis) + + +###################################################################### +# Overlaying Image Axes +# --------------------- +# +# For the common use case of plotting image axis on an ``ArrayDisplay``, +# the ``set_line_hillas()`` method is provided for convenience. The +# following example shows its use: +# + + +input_url = "dataset://gamma_LaPalma_baseline_20Zd_180Az_prod3b_test.simtel.gz" + + +###################################################################### +# First, we define a function to plot the array with overlaid lines for +# the image axes +# + + +def plot_event(event, subarray, ax): + """ + Draw an ArrayDisplay with image axes and the + true and reconstructed impact position overlaid + """ + + event.pointing.array_azimuth + disp = ArrayDisplay(subarray, axes=ax) + + hillas_dict = {tid: tel.parameters.hillas for tid, tel in event.dl1.tel.items()} + core_dict = {tid: tel.parameters.core.psi for tid, tel in event.dl1.tel.items()} + + disp.set_line_hillas( + hillas_dict, + core_dict, + 500, + ) + + reco_shower = event.dl2.stereo.geometry["HillasReconstructor"] + + ax.scatter( + event.simulation.shower.core_x, + event.simulation.shower.core_y, + s=200, + c="k", + marker="x", + label="True Impact", + ) + ax.scatter( + reco_shower.core_x, + reco_shower.core_y, + s=200, + c="r", + marker="x", + label="Estimated Impact", + ) + + ax.legend() + + +###################################################################### +# Now, we can loop through some events and plot them. Here we apply +# default calibration, image processing, and reconstruction, however it is +# better to use ``ctapipe-process`` with a well-defined configuration to +# do this in reality. Note that some events will not have images bright +# enough to do parameterization or reconstruction, so they will have no +# image axis lines or no estimated impact position. +# + +fig, ax = plt.subplots(5, 3, figsize=(20, 40), constrained_layout=True) +ax = ax.ravel() + +with EventSource(input_url, max_events=15, focal_length_choice="EQUIVALENT") as source: + calib = CameraCalibrator(subarray=source.subarray) + process_images = ImageProcessor(subarray=source.subarray) + process_shower = ShowerProcessor(subarray=source.subarray) + + for i, event in enumerate(source): + calib(event) + process_images(event) + process_shower(event) + plot_event(event, source.subarray, ax=ax[i]) diff --git a/examples/visualization/camera_display.py b/examples/visualization/camera_display.py new file mode 100644 index 00000000000..01d7a18ee36 --- /dev/null +++ b/examples/visualization/camera_display.py @@ -0,0 +1,410 @@ +""" +Displaying Camera Images +======================== + +""" + +import astropy.coordinates as c +import astropy.units as u +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.animation import FuncAnimation +from matplotlib.colors import PowerNorm + +from ctapipe.coordinates import CameraFrame, EngineeringCameraFrame, TelescopeFrame +from ctapipe.image import hillas_parameters, tailcuts_clean, toymodel +from ctapipe.instrument import SubarrayDescription +from ctapipe.io import EventSource +from ctapipe.visualization import CameraDisplay + +###################################################################### +# First, let’s create a fake Cherenkov image from a given +# ``CameraGeometry`` and fill it with some data that we can draw later. +# + +# load an example camera geometry from a simulation file +subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst") +geom = subarray.tel[100].camera.geometry + +# create a fake camera image to display: +model = toymodel.Gaussian( + x=0.2 * u.m, + y=0.0 * u.m, + width=0.05 * u.m, + length=0.15 * u.m, + psi="35d", +) + +image, sig, bg = model.generate_image(geom, intensity=1500, nsb_level_pe=10) +mask = tailcuts_clean(geom, image, picture_thresh=15, boundary_thresh=5) + +###################################################################### +geom + + +###################################################################### +# Displaying Images +# ----------------- +# +# The simplest plot is just to generate a CameraDisplay with an image in +# its constructor. A figure and axis will be created automatically +# + +CameraDisplay(geom) + + +###################################################################### +# You can also specify the initial ``image``, ``cmap`` and ``norm`` +# (colomap and normalization, see below), ``title`` to use. You can +# specify ``ax`` if you want to draw the camera on an existing +# *matplotlib* ``Axes`` object (otherwise one is created). +# +# To change other options, or to change options dynamically, you can call +# the relevant functions of the ``CameraDisplay`` object that is returned. +# For example to add a color bar, call ``add_colorbar()``, or to change +# the color scale, modify the ``cmap`` or ``norm`` properties directly. +# + + +###################################################################### +# Choosing a coordinate frame +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# The ``CameraGeometry`` object contains a ``ctapipe.coordinates.Frame`` +# used by ``CameraDisplay`` to draw the camera in the correct orientation +# and distance units. The default frame is the ``CameraFrame``, which will +# display the camera in units of *meters* and with an orientation that the +# top of the camera (when parked) is aligned to the X-axis. To show the +# camera in another orientation, it’s useful to apply a coordinate +# transform to the ``CameraGeometry`` before passing it to the +# ``CameraDisplay``. The following ``Frames`` are supported: \* +# ``EngineeringCameraFrame`` : similar to CameraFrame, but with the top of +# the camera aligned to the Y axis \* ``TelescopeFrame``: In *degrees* (on +# the sky) coordinates relative to the telescope Alt/Az pointing position, +# with the Alt axis pointing upward. +# + +fig, ax = plt.subplots(1, 3, figsize=(15, 4)) +CameraDisplay(geom, image=image, ax=ax[0]) +CameraDisplay(geom.transform_to(EngineeringCameraFrame()), image=image, ax=ax[1]) +CameraDisplay(geom.transform_to(TelescopeFrame()), image=image, ax=ax[2]) + + +###################################################################### +# Note the the name of the Frame appears in the lower-right corner +# + + +###################################################################### +# For the rest of this demo, let’s use the ``TelescopeFrame`` +# + +geom_camframe = geom +geom = geom_camframe.transform_to(EngineeringCameraFrame()) + + +###################################################################### +# Changing the color map and scale +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# + + +###################################################################### +# CameraDisplay supports any `matplotlib color +# map `__ +# It is **highly recommended** to use a *perceptually uniform* map, unless +# you have a good reason not to. +# + +fig, ax = plt.subplots(1, 3, figsize=(15, 4)) +for ii, cmap in enumerate(["PuOr_r", "rainbow", "twilight"]): + disp = CameraDisplay(geom, image=image, ax=ax[ii], title=cmap) + disp.add_colorbar() + disp.cmap = cmap + + +###################################################################### +# By default the minimum and maximum of the color bar are set +# automatically by the data in the image. To choose fixed limits, use:\` +# + +fig, ax = plt.subplots(1, 3, figsize=(15, 4)) +for ii, minmax in enumerate([(10, 50), (-10, 10), (1, 100)]): + disp = CameraDisplay(geom, image=image, ax=ax[ii], title=minmax) + disp.add_colorbar() + disp.set_limits_minmax(minmax[0], minmax[1]) + + +###################################################################### +# Or you can set the maximum limit by percentile of the charge +# distribution: +# + +fig, ax = plt.subplots(1, 3, figsize=(15, 4)) +for ii, pct in enumerate([30, 50, 90]): + disp = CameraDisplay(geom, image=image, ax=ax[ii], title=f"{pct} %") + disp.add_colorbar() + disp.set_limits_percent(pct) + + +###################################################################### +# Using different normalizations +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# + + +###################################################################### +# You can choose from several preset normalizations (lin, log, symlog) and +# also provide a custom normalization, for example a ``PowerNorm``: +# + + +fig, axes = plt.subplots(2, 2, figsize=(14, 10)) +norms = ["lin", "log", "symlog", PowerNorm(0.5)] + +for norm, ax in zip(norms, axes.flatten()): + disp = CameraDisplay(geom, image=image, ax=ax) + disp.norm = norm + disp.add_colorbar() + ax.set_title(str(norm)) + +axes[1, 1].set_title("PowerNorm(0.5)") +plt.show() + + +###################################################################### +# Overlays +# -------- +# + + +###################################################################### +# Marking pixels +# ~~~~~~~~~~~~~~ +# +# here we will mark pixels in the image mask. That will change their +# outline color +# + +fig, ax = plt.subplots(1, 2, figsize=(10, 4)) +disp = CameraDisplay( + geom, image=image, cmap="gray", ax=ax[0], title="Image mask in green" +) +disp.highlight_pixels(mask, alpha=0.8, linewidth=2, color="green") + +disp = CameraDisplay( + geom, image=image, cmap="gray", ax=ax[1], title="Image mask in green (zoom)" +) +disp.highlight_pixels(mask, alpha=1, linewidth=3, color="green") + +ax[1].set_ylim(-0.5, 0.5) +ax[1].set_xlim(-0.5, 0.5) + + +###################################################################### +# Drawing a Hillas-parameter ellipse +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# + + +###################################################################### +# For this, we will first compute some Hillas Parameters in the current +# frame: +# + +clean_image = image.copy() +clean_image[~mask] = 0 +hillas = hillas_parameters(geom, clean_image) + +plt.figure(figsize=(6, 6)) +disp = CameraDisplay(geom, image=image, cmap="gray_r") +disp.highlight_pixels(mask, alpha=0.5, color="dodgerblue") +disp.overlay_moments(hillas, color="red", linewidth=3, with_label=False) + + +###################################################################### +# Drawing a marker at a coordinate +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# This depends on the coordinate frame of the ``CameraGeometry``. Here we +# will sepcify the coordinate the ``EngineerngCameraFrame``, but if you +# have enough information to do the coordinate transform, you could use +# ``ICRS`` coordinates and overlay star positions. ``CameraDisplay`` will +# convert the coordinate you pass in to the ``Frame`` of the display +# automatically (if sufficient frame attributes are set). +# +# Note that the parameter ``keep_old`` is False by default, meaning adding +# a new point will clear the previous ones (useful for animations, but +# perhaps unexpected for a static plot). Set it to ``True`` to plot +# multiple markers. +# + +plt.figure(figsize=(6, 6)) +disp = CameraDisplay(geom, image=image, cmap="gray_r") + +coord = c.SkyCoord(x=0.5 * u.m, y=0.7 * u.m, frame=geom.frame) +coord_in_another_frame = c.SkyCoord(x=0.5 * u.m, y=0.7 * u.m, frame=CameraFrame()) +disp.overlay_coordinate(coord, markersize=20, marker="*") +disp.overlay_coordinate( + coord_in_another_frame, markersize=20, marker="*", keep_old=True +) + + +###################################################################### +# Generating an animation +# ----------------------- +# +# Here we will make an animation of fake events by re-using a single +# display (much faster than generating a new one each time) +# + + +subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst") +geom = subarray.tel[1].camera.geometry + +fov = 1.0 +maxwid = 0.05 +maxlen = 0.1 + +fig, ax = plt.subplots(1, 1, figsize=(8, 6)) +disp = CameraDisplay(geom, ax=ax) # we only need one display (it can be re-used) +disp.cmap = "inferno" +disp.add_colorbar(ax=ax) + + +def update(frame): + """this function will be called for each frame of the animation""" + x, y = np.random.uniform(-fov, fov, size=2) + width = np.random.uniform(0.01, maxwid) + length = np.random.uniform(width, maxlen) + angle = np.random.uniform(0, 180) + intens = width * length * (5e4 + 1e5 * np.random.exponential(2)) + + model = toymodel.Gaussian( + x=x * u.m, + y=y * u.m, + width=width * u.m, + length=length * u.m, + psi=angle * u.deg, + ) + image, _, _ = model.generate_image( + geom, + intensity=intens, + nsb_level_pe=5, + ) + disp.image = image + + +# Create the animation and convert to a displayable video: +anim = FuncAnimation(fig, func=update, frames=10, interval=200) +plt.show() + +###################################################################### +# Using CameraDisplays interactively +# ---------------------------------- +# +# ``CameraDisplays`` can be used interactivly whe displayed in a window, +# and also when using Jupyter notebooks/lab with appropriate backends. + +###################################################################### +# When this is the case, the same ``CameraDisplay`` object can be re-used. +# We can’t show this here in the documentation, but creating an animation +# when in a matplotlib window is quite easy! Try this in an interactive +# ipython session: +# + +###################################################################### +# Running interactive displays in a matplotlib window +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# ipython -i --maplotlib=auto + + +###################################################################### +# That will open an ipython session with matplotlib graphics in a separate +# thread, meaning that you can type code and interact with plots +# simultaneneously. +# +# In the ipython session try running the following code and you will see +# an animation (here in the documentation, it will of course be static) +# +# First we load some real data so we have a nice image to view: +# + + +###################################################################### +DATA = "dataset://gamma_20deg_0deg_run1___cta-prod5-lapalma_desert-2158m-LaPalma-dark_100evts.simtel.zst" + +with EventSource( + DATA, + max_events=1, + focal_length_choice="EQUIVALENT", +) as source: + event = next(iter(source)) + +tel_id = list(event.r0.tel.keys())[0] +geom = source.subarray.tel[tel_id].camera.geometry +waveform = event.r0.tel[tel_id].waveform +n_chan, n_pix, n_samp = waveform.shape + + +###################################################################### +# Running the following the will bring up a window and animate the shower +# image as a function of time. +# + +disp = CameraDisplay(geom) + +for ii in range(n_samp): + disp.image = waveform[0, :, ii] + plt.pause(0.1) # this lets matplotlib re-draw the scene + + +###################################################################### +# The output will be similar to the static animation created as follows: +# + +fig, ax = plt.subplots(1, 1) +disp = CameraDisplay(geom, ax=ax) +disp.add_colorbar() +disp.autoscale = False + + +def draw_sample(frame): + ax.set_title(f"sample: {frame}") + disp.set_limits_minmax(200, 400) + disp.image = waveform[0, :, frame] + + +anim = FuncAnimation(fig, func=draw_sample, frames=n_samp, interval=100) +plt.show() + +###################################################################### +# Making it clickable +# ~~~~~~~~~~~~~~~~~~~ +# +# Also when running in a window, you can enable the +# ``disp.enable_pixel_picker()`` option. This will then allow the user to +# click a pixel and a function will run. By default the function simply +# prints the pixel and value to stdout, however you can override the +# function ``on_pixel_clicked(pix_id)`` to do anything you want by making +# a subclass +# + + +class MyCameraDisplay(CameraDisplay): + def on_pixel_clicked(self, pix_id): + print(f"{pix_id=} has value {self.image[pix_id]:.2f}") + + +disp = MyCameraDisplay(geom, image=image) +disp.enable_pixel_picker() + + +###################################################################### +# then, when a user clicks a pixel it would print: +# +# :: +# +# pixel 5 has value 2.44 +# diff --git a/pyproject.toml b/pyproject.toml index 7501b7da070..73a42e69d99 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,6 +28,11 @@ filterwarnings = [ "error::astropy.utils.exceptions.AstropyDeprecationWarning", "error::ctapipe.utils.deprecation.CTAPipeDeprecationWarning", ] +norecursedirs = [ + ".git", + "_build", + "auto_examples", +] [tool.towncrier] package = "ctapipe" diff --git a/setup.cfg b/setup.cfg index 767a6b458f6..dbd35656b60 100644 --- a/setup.cfg +++ b/setup.cfg @@ -61,12 +61,14 @@ docs = nbsphinx numpydoc sphinx-design + sphinx-gallery jupyter notebook graphviz pandas ipython ffmpeg-python + pypandoc dev = setuptools_scm[toml]