Array plotting

ctapipe/image/tests/test_timing_parameters.py

@@ -13,8 +13,8 @@ def test_grad_fit():
     intercept = 1
 
     timing = timing_parameters(
-        pix_x=np.zeros(4) * u.deg,
-        pix_y=np.arange(4) * u.deg,
+        pix_x=np.arange(4) * u.deg,
+        pix_y=np.zeros(4) * u.deg,
         image=np.ones(4),
         peak_time=intercept * u.ns + grad * np.arange(4) * u.ns,
         rotation_angle=0 * u.deg
@@ -27,8 +27,8 @@ def test_grad_fit():
     # Then try a different rotation angle
     rot_angle = 20 * u.deg
     timing_rot20 = timing_parameters(
-        pix_x=np.zeros(4) * u.deg,
-        pix_y=np.arange(4) * u.deg,
+        pix_x=np.arange(4) * u.deg,
+        pix_y=np.zeros(4) * u.deg,
         image=np.ones(4),
         peak_time=intercept * u.ns +
         grad * np.arange(4) * u.ns,

ctapipe/image/timing_parameters.py

@@ -23,7 +23,7 @@ class TimingParameterizationError(RuntimeError):
     pass
 
 
-def rotate_translate(pixel_pos_x, pixel_pos_y, phi):
+def rotate_translate(pixel_pos_x, pixel_pos_y, psi):
     """
     Function to perform rotation and translation of pixel lists
 
@@ -41,9 +41,10 @@ def rotate_translate(pixel_pos_x, pixel_pos_y, phi):
         ndarray,ndarray: Transformed pixel x and y coordinates
 
     """
-
-    pixel_pos_rot_x = pixel_pos_x * np.cos(phi) - pixel_pos_y * np.sin(phi)
-    pixel_pos_rot_y = pixel_pos_x * np.sin(phi) + pixel_pos_y * np.cos(phi)
+    # Here we are derotating, so we use a minus sign
+    psi = -psi
+    pixel_pos_rot_x = pixel_pos_x * np.cos(psi) - pixel_pos_y * np.sin(psi)
+    pixel_pos_rot_y = pixel_pos_x * np.sin(psi) + pixel_pos_y * np.cos(psi)
     return pixel_pos_rot_x, pixel_pos_rot_y
 
 
@@ -62,7 +63,7 @@ def timing_parameters(pix_x, pix_y, image, peak_time, rotation_angle):
     peak_time : array_like
         Pixel times corresponding
     rotation_angle: float
-        Rotation angle fo the image major axis
+        Rotation angle for the image major axis, namely psi (in radians)
 
     Returns
     -------
@@ -75,13 +76,31 @@ def timing_parameters(pix_x, pix_y, image, peak_time, rotation_angle):
     image = np.asanyarray(image, dtype=np.float64)
     peak_time = np.asanyarray(peak_time, dtype=np.float64)
 
+    # select only the pixels in the cleaned image that are greater than zero.
+    # This is to allow to use a dilated mask (which might be better):
+    # we need to exclude possible pixels with zero signal after cleaning.
+
+    mask = np.ma.masked_where(image > 0, image).mask
+    pix_x = pix_x[mask]
+    pix_y = pix_y[mask]
+    image = image[mask]
+    peak_time = peak_time[mask]
+
     assert pix_x.shape == image.shape
     assert pix_y.shape == image.shape
     assert peak_time.shape == image.shape
 
-    # Rotate pixels by our image axis
-    pix_x_rot, pix_y_rot = rotate_translate(pix_x, pix_y, rotation_angle)
-    gradient, intercept = np.polyfit(pix_y_rot, peak_time, deg=1, w=np.sqrt(image))
+    # since the values of peak_pos are integers, sometimes the mask is constant
+    # for all the selected pixels. Asking for a mask to have at least 3 different
+    # values for the peak_pos in the selected pixels seems reasonable.
+
+    if np.unique(peak_time).size > 2:
+        pix_x_rot, pix_y_rot = rotate_translate(pix_x, pix_y, rotation_angle)
+        gradient, intercept = np.polyfit(x=pix_x_rot, y=peak_time,
+                                         deg=1, w=np.sqrt(image))
+    else:
+        gradient = 0.
+        intercept = 0.
 
     return TimingParameters(gradient=gradient * (peak_time.unit / unit),
                             intercept=intercept * unit)

ctapipe/visualization/mpl_array.py

@@ -138,6 +138,8 @@ def set_vector_uv(self, u, v, c=None, **kwargs):
             x-component of direction vector
         v: array[num_tels]
             y-component of direction vector
+        c: color or list of colors
+            vector color for each telescope (or one for all)
         kwargs:
             extra args passed to plt.quiver(), ignored on subsequent updates
         """
@@ -150,6 +152,9 @@ def set_vector_uv(self, u, v, c=None, **kwargs):
                 coords.x, coords.y,
                 u, v,
                 color=c,
+                scale_units='xy',
+                angles='xy',
+                scale=1,
                 **kwargs
             )
         else:
@@ -171,26 +176,47 @@ def set_vector_rho_phi(self, rho, phi, c=None, **kwargs):
         u, v = polar_to_cart(rho, phi)
         self.set_vector_uv(u, v, c=c, **kwargs)
 
-    def set_vector_hillas(self, hillas_dict, angle_offset=180 * u.deg):
+    def set_vector_hillas(self, hillas_dict, length, time_gradient, angle_offset):
         """
-        helper function to set the vector angle and length from a set of
-        Hillas parameters.
+        Function to set the vector angle and length from a set of Hillas parameters.
+
+        In order to proper use the arrow on the ground, also a dictionary with the time
+        gradients for the different telescopes is needed. If the gradient is 0 the arrow
+        is not plotted on the ground, whereas if the value of the gradient is negative,
+        the arrow is rotated by 180 degrees (Angle(angle_offset) not added).
+
+        This plotting behaviour has been tested with the timing_parameters function
+        in ctapipe/image.
 
         Parameters
         ----------
         hillas_dict: Dict[int, HillasParametersContainer]
             mapping of tel_id to Hillas parameters
+        length: Float
+            length of the arrow (in meters)
+        time_gradient: Dict[int, value of time gradient (no units)]
+            dictionary for value of the time gradient for each telescope
+        angle_offset: Float
+            This should be the event.mcheader.run_array_direction[0] parameter
+
         """
 
         # rot_angle_ellipse is psi parameter in HillasParametersContainer
-
         rho = np.zeros(self.subarray.num_tels) * u.m
         rot_angle_ellipse = np.zeros(self.subarray.num_tels) * u.deg
 
         for tel_id, params in hillas_dict.items():
             idx = self.subarray.tel_indices[tel_id]
-            rho[idx] = 1.0 * u.m  # params.length
-            rot_angle_ellipse[idx] = Angle(params.psi)+ Angle(angle_offset)
+            rho[idx] = length * u.m
+
+            if time_gradient[tel_id] > 0.01:
+                params.psi = Angle(params.psi)
+                angle_offset = Angle(angle_offset)
+                rot_angle_ellipse[idx] = params.psi + angle_offset + 180 * u.deg
+            elif time_gradient[tel_id] < -0.01:
+                rot_angle_ellipse[idx] = params.psi + angle_offset
+            else:
+                rho[idx] = 0 * u.m
 
         self.set_vector_rho_phi(rho=rho, phi=rot_angle_ellipse)
 
@@ -217,7 +243,7 @@ def set_line_hillas(self, hillas_dict, range, **kwargs):
             y_0 = coords[idx].y.value
             m = np.tan(Angle(params.psi))
             x = x_0 + np.linspace(-range, range, 50)
-            y = y_0 + m * (x-x_0)
+            y = y_0 + m * (x - x_0)
             distance = np.sqrt((x - x_0) ** 2 + (y - y_0) ** 2)
             mask = np.ma.masked_where(distance < range, distance).mask
             self.axes.plot(x[mask], y[mask], color=c[idx], **kwargs)
@@ -228,7 +254,7 @@ def add_labels(self):
         py = self.tel_coords.y.value
         for tel, x, y in zip(self.subarray.tels, px, py):
             name = str(tel)
-            lab = self.axes.text(x, y, name, fontsize=8)
+            lab = self.axes.text(x, y, name, fontsize=8, clip_on=True)
             self._labels.append(lab)
 
     def remove_labels(self):

ctapipe/visualization/tests/test_mpl.py

@@ -6,7 +6,7 @@
 from ctapipe.instrument import (CameraGeometry, SubarrayDescription,
                                 TelescopeDescription)
 from ctapipe.io.containers import HillasParametersContainer
-from numpy import ones
+from numpy import ones, zeros, arange
 from astropy import units as u
 
 
@@ -55,6 +55,7 @@ def test_camera_display_multiple():
 
 def test_array_display():
     from ctapipe.visualization.mpl_array import ArrayDisplay
+    from ctapipe.image.timing_parameters import timing_parameters
 
     # build a test subarray:
     tels = dict()
@@ -80,10 +81,30 @@ def test_array_display():
 
     # test using hillas params:
     hillas_dict = {
-        1: HillasParametersContainer(length=1.0 * u.m, psi=90 * u.deg),
-        2: HillasParametersContainer(length=200 * u.cm, psi="95deg"),
+        1: HillasParametersContainer(length=100.0 * u.m, psi=90 * u.deg),
+        2: HillasParametersContainer(length=20000 * u.cm, psi="95deg"),
     }
-    ad.set_vector_hillas(hillas_dict)
+
+    grad = 2
+    intercept = 1
+
+    rot_angle = 20 * u.deg
+    timing_rot20 = timing_parameters(
+        pix_x=arange(4) * u.deg,
+        pix_y=zeros(4) * u.deg,
+        image=ones(4),
+        peak_time=intercept * u.ns + grad * arange(4) * u.ns,
+        rotation_angle=rot_angle
+    )
+    gradient_dict = {
+        1: timing_rot20.gradient.value,
+        2: timing_rot20.gradient.value,
+    }
+    ad.set_vector_hillas(hillas_dict=hillas_dict,
+                         length=500,
+                         time_gradient=gradient_dict,
+                         angle_offset=0 * u.deg)
+
     ad.set_line_hillas(hillas_dict, range=300)
     ad.add_labels()
     ad.remove_labels()