DL1 writer tool (ctapipe-stage1-process)

ci/install.sh

@@ -23,9 +23,9 @@ if [[ "$CONDA" == "true" ]]; then
     # Useful for debugging any issues with conda
     conda info -a
 
-    sed -i -e "s/- python=.*/- python=$PYTHON_VERSION/g" environment.yml
-    travis_wait 20 conda env create -n cta-dev --file environment.yml
-    conda activate cta-dev
+	sed -i -e "s/- python=.*/- python=$TRAVIS_PYTHON_VERSION/g" environment.yml
+	travis_wait 20 conda env create -n cta-dev --file environment.yml
+	conda activate cta-dev
 else
     pip install -U pip
 fi

ctapipe/containers.py

@@ -83,20 +83,19 @@ class EventType(enum.Enum):
 
 class EventIndexContainer(Container):
     """ index columns to include in event lists, common to all data levels"""
-
     container_prefix = ""  # don't want to prefix these
-
-    event_id = Field(0, "event identifier")
     obs_id = Field(0, "observation identifier")
+    event_id = Field(0, "event identifier")
 
 
-class TelEventIndexContainer(EventIndexContainer):
+class TelEventIndexContainer(Container):
     """
     index columns to include in telescope-wise event lists, common to all data
     levels that have telescope-wise information
     """
-
     container_prefix = ""  # don't want to prefix these
+    obs_id = Field(0, "observation identifier")
+    event_id = Field(0, "event identifier")
     tel_id = Field(0, "telescope identifier")
 
 
@@ -181,11 +180,30 @@ class TimingParametersContainer(Container):
 class MorphologyContainer(Container):
     """ Parameters related to pixels surviving image cleaning """
 
-    num_pixels = Field(nan, "Number of usable pixels")
-    num_islands = Field(nan, "Number of distinct islands in the image")
-    num_small_islands = Field(nan, "Number of <= 2 pixel islands")
-    num_medium_islands = Field(nan, "Number of 2-50 pixel islands")
-    num_large_islands = Field(nan, "Number of > 10 pixel islands")
+    num_pixels = Field(-1, "Number of usable pixels")
+    num_islands = Field(-1, "Number of distinct islands in the image")
+    num_small_islands = Field(-1, "Number of <= 2 pixel islands")
+    num_medium_islands = Field(-1, "Number of 2-50 pixel islands")
+    num_large_islands = Field(-1, "Number of > 50 pixel islands")
+
+
+class StatisticsContainer(Container):
+    """Store descriptive statistics"""
+
+    max = Field(nan, "value of pixel with maximum intensity")
+    min = Field(nan, "value of pixel with minimum intensity")
+    mean = Field(nan, "mean intensity")
+    std = Field(nan, "standard deviation of intensity")
+    skewness = Field(nan, "skewness of intensity")
+    kurtosis = Field(nan, "kurtosis of intensity")
+
+
+class IntensityStatisticsContainer(StatisticsContainer):
+    container_prefix = "intensity"
+
+
+class PeakTimeStatisticsContainer(StatisticsContainer):
+    container_prefix = "peak_time"
 
 
 class ImageParametersContainer(Container):
@@ -197,6 +215,12 @@ class ImageParametersContainer(Container):
     leakage = Field(LeakageContainer(), "Leakage Parameters")
     concentration = Field(ConcentrationContainer(), "Concentration Parameters")
     morphology = Field(MorphologyContainer(), "Image Morphology Parameters")
+    intensity_statistics = Field(
+        IntensityStatisticsContainer(), "Intensity image statistics"
+    )
+    peak_time_statistics = Field(
+        PeakTimeStatisticsContainer(), "Peak time image statistics"
+    )
 
 
 class DL1CameraContainer(Container):
@@ -230,7 +254,7 @@ class DL1CameraContainer(Container):
     parameters = Field(ImageParametersContainer(), "Parameters derived from images")
 
 
-class MCDL1CameraContainer(DL1CameraContainer):
+class MCDL1CameraContainer(Container):
     """ Contains all fields of the DL1CameraContainer, but adds fields for simulated
     DL1 image information."""
 

ctapipe/core/tool.py

@@ -214,7 +214,7 @@ def run(self, argv=None):
             Provenance().start_activity(self.name)
             self.setup()
             self.is_setup = True
-            self.log.info(f"CONFIG: {self.get_current_config()}")
+            self.log.debug(f"CONFIG: {self.get_current_config()}")
             Provenance().add_config(self.get_current_config())
             self.start()
             self.finish()

ctapipe/image/__init__.py

@@ -1,9 +1,20 @@
-from .hillas import *
+from .hillas import (
+    hillas_parameters, HillasParameterizationError, camera_to_shower_coordinates,
+)
+from .timing import timing_parameters
+from .leakage import leakage
+from .concentration import concentration
+from .statistics import descriptive_statistics
+from .morphology import (
+    number_of_islands,
+    number_of_island_sizes,
+    morphology_parameters,
+    largest_island,
+)
+
 from .cleaning import *
 from .pixel_likelihood import *
 from .extractor import *
 from .reducer import *
-from .muon import *
 from .geometry_converter import *
-from .leakage import *
-from .concentration import concentration
+from .muon import *

ctapipe/image/cleaning.py

@@ -8,15 +8,13 @@
     "mars_cleaning_1st_pass",
     "fact_image_cleaning",
     "apply_time_delta_cleaning",
-    "number_of_islands",
     "ImageCleaner",
     "TailcutsImageCleaner",
 ]
 
 from abc import abstractmethod
 
 import numpy as np
-from scipy.sparse.csgraph import connected_components
 
 from ..core.component import TelescopeComponent
 from ..core.traits import (
@@ -202,42 +200,6 @@ def dilate(geom, mask):
     return mask | geom.neighbor_matrix_sparse.dot(mask)
 
 
-def number_of_islands(geom, mask):
-    """
-    Search a given pixel mask for connected clusters.
-    This can be used to seperate between gamma and hadronic showers.
-
-    Parameters
-    ----------
-    geom: `~ctapipe.instrument.CameraGeometry`
-        Camera geometry information
-    mask: ndarray
-        input mask (array of booleans)
-
-    Returns
-    -------
-    num_islands: int
-        Total number of clusters
-    island_labels: ndarray
-        Contains cluster membership of each pixel.
-        Dimension equals input geometry.
-        Entries range from 0 (not in the pixel mask) to num_islands.
-    """
-    # compress sparse neighbor matrix
-    neighbor_matrix_compressed = geom.neighbor_matrix_sparse[mask][:, mask]
-    # pixels in no cluster have label == 0
-    island_labels = np.zeros(geom.n_pixels, dtype="int32")
-
-    num_islands, island_labels_compressed = connected_components(
-        neighbor_matrix_compressed, directed=False
-    )
-
-    # count clusters from 1 onwards
-    island_labels[mask] = island_labels_compressed + 1
-
-    return num_islands, island_labels
-
-
 def apply_time_delta_cleaning(
     geom, mask, arrival_times, min_number_neighbors, time_limit
 ):
@@ -367,30 +329,6 @@ def fact_image_cleaning(
     return pixels_to_keep
 
 
-def largest_island(islands_labels):
-    """Find the biggest island and filter it from the image.
-
-    This function takes a list of islands in an image and isolates the largest one
-    for later parametrization.
-
-    Parameters
-    ----------
-    islands_labels : array
-        Flattened array containing a list of labelled islands from a cleaned image.
-        Second value returned by the function 'number_of_islands'.
-
-    Returns
-    -------
-    islands_labels : array
-        A boolean mask created from the input labels and filtered for the largest island.
-        If no islands survived the cleaning the array is all False.
-
-    """
-    if np.count_nonzero(islands_labels) == 0:
-        return np.zeros(islands_labels.shape, dtype="bool")
-    return islands_labels == np.argmax(np.bincount(islands_labels[islands_labels > 0]))
-
-
 class ImageCleaner(TelescopeComponent):
     """
     Abstract class for all configurable Image Cleaning algorithms.   Use

ctapipe/image/extractor.py

@@ -26,8 +26,8 @@
 from ctapipe.core import TelescopeComponent
 from numba import njit, prange, guvectorize, float64, float32, int64
 
-from .cleaning import number_of_islands, largest_island, tailcuts_clean
-from .timing_parameters import timing_parameters
+from . import number_of_islands, largest_island, tailcuts_clean
+from .timing import timing_parameters
 from .hillas import hillas_parameters, camera_to_shower_coordinates
 
 

ctapipe/image/morphology.py

@@ -0,0 +1,126 @@
+import numpy as np
+from scipy.sparse.csgraph import connected_components
+
+from ..containers import MorphologyContainer
+
+
+def number_of_islands(geom, mask):
+    """
+    Search a given pixel mask for connected clusters.
+    This can be used to seperate between gamma and hadronic showers.
+
+    Parameters
+    ----------
+    geom: `~ctapipe.instrument.CameraGeometry`
+        Camera geometry information
+    mask: ndarray
+        input mask (array of booleans)
+
+    Returns
+    -------
+    num_islands: int
+        Total number of clusters
+    island_labels: ndarray
+        Contains cluster membership of each pixel.
+        Dimension equals input geometry.
+        Entries range from 0 (not in the pixel mask) to num_islands.
+    """
+    # compress sparse neighbor matrix
+    neighbor_matrix_compressed = geom.neighbor_matrix_sparse[mask][:, mask]
+    # pixels in no cluster have label == 0
+    island_labels = np.zeros(geom.n_pixels, dtype="int32")
+
+    num_islands, island_labels_compressed = connected_components(
+        neighbor_matrix_compressed, directed=False
+    )
+
+    # count clusters from 1 onwards
+    island_labels[mask] = island_labels_compressed + 1
+
+    return num_islands, island_labels
+
+
+def number_of_island_sizes(island_labels):
+    '''
+    Return number of small, medium and large islands
+
+    Parameters
+    ----------
+    island_labels: array[int]
+        Array with island labels, (second return value of ``number_of_islands``)
+
+    Returns
+    -------
+    n_small: int
+        number of islands with less than 3 pixels
+    n_medium: int
+        number of islands with 3 <= n_pixels <= 50
+    n_large: int
+        number of islands with more than 50 pixels
+    '''
+
+    # count number of pixels in each island, remove 0 = no island
+    island_sizes = np.bincount(island_labels)[1:]
+
+    # remove islands of size 0 (if labels are not consecutive)
+    # should not happen, but easy to check
+    island_sizes = island_sizes[island_sizes > 0]
+
+    small = island_sizes <= 2
+    large = island_sizes > 50
+    n_medium = np.count_nonzero(~(small | large))
+
+    return np.count_nonzero(small), n_medium, np.count_nonzero(large)
+
+
+def largest_island(islands_labels):
+    """Find the biggest island and filter it from the image.
+
+    This function takes a list of islands in an image and isolates the largest one
+    for later parametrization.
+
+    Parameters
+    ----------
+    islands_labels : array
+        Flattened array containing a list of labelled islands from a cleaned image.
+        Second value returned by the function 'number_of_islands'.
+
+    Returns
+    -------
+    islands_labels : array
+        A boolean mask created from the input labels and filtered for the largest island.
+        If no islands survived the cleaning the array is all False.
+
+    """
+    if np.count_nonzero(islands_labels) == 0:
+        return np.zeros(islands_labels.shape, dtype="bool")
+    return islands_labels == np.argmax(np.bincount(islands_labels[islands_labels > 0]))
+
+
+def morphology_parameters(geom, image_mask) -> MorphologyContainer:
+    """
+    Compute image morphology parameters
+
+    Parameters
+    ----------
+    geom: ctapipe.instrument.camera.CameraGeometry
+        camera description
+    image_mask: np.ndarray(bool)
+       image of pixels surviving cleaning (True=survives)
+
+    Returns
+    -------
+    MorphologyContainer: parameters related to the morphology
+    """
+
+    num_islands, island_labels = number_of_islands(geom=geom, mask=image_mask)
+
+    n_small, n_medium, n_large = number_of_island_sizes(island_labels)
+
+    return MorphologyContainer(
+        num_pixels=np.count_nonzero(image_mask),
+        num_islands=num_islands,
+        num_small_islands=n_small,
+        num_medium_islands=n_medium,
+        num_large_islands=n_large,
+    )

ctapipe/image/statistics.py

@@ -0,0 +1,15 @@
+from scipy.stats import skew, kurtosis
+
+from ..containers import StatisticsContainer
+
+
+def descriptive_statistics(values, container_class=StatisticsContainer) -> StatisticsContainer:
+    """ compute intensity statistics of an image  """
+    return container_class(
+        max=values.max(),
+        min=values.min(),
+        mean=values.mean(),
+        std=values.std(),
+        skewness=skew(values),
+        kurtosis=kurtosis(values),
+    )