Repository: cta-observatory/ctapipe (python)

Issue:
CameraGeometry position_to_pix_index returns -1 for nan input
https://github.com/cta-observatory/ctapipe/blob/39762693534f7dc9b884684c3f82e153bd9da560/ctapipe/instrument/camera/geometry.py#L859

For example, of one uses it on image centroid coordinates to make a "pixelized" 2d histogram, all non-reconstructed images end up in the highest-index pixel (=-1). It would be better to return np.nan.

Pull request description:
Use intmin for invalid pixel positions, allow empty arguments
Fixes #1828

Patch:
diff --git a/ctapipe/instrument/camera/geometry.py b/ctapipe/instrument/camera/geometry.py
index dd55b144..1d83a1ba 100644
--- a/ctapipe/instrument/camera/geometry.py
+++ b/ctapipe/instrument/camera/geometry.py
@@ -954,12 +954,13 @@ class CameraGeometry:
         pix_indices: Pixel index or array of pixel indices. Returns -1 if position falls
                     outside camera
         """
-
         if not self._all_pixel_areas_equal:
             logger.warning(
                 " Method not implemented for cameras with varying pixel sizes"
             )
         unit = x.unit
+        scalar = x.ndim == 0
+
         points_searched = np.dstack([x.to_value(unit), y.to_value(unit)])
         circum_rad = self._pixel_circumradius[0].to_value(unit)
         kdtree = self._kdtree
@@ -969,8 +970,9 @@ class CameraGeometry:
         del dist
         pix_indices = pix_indices.flatten()
 
+        invalid = np.iinfo(pix_indices.dtype).min
         # 1. Mark all points outside pixel circumeference as lying outside camera
-        pix_indices[pix_indices == self.n_pixels] = -1
+        pix_indices[pix_indices == self.n_pixels] = invalid
 
         # 2. Accurate check for the remaing cases (within circumference, but still outside
         # camera). It is first checked if any border pixel numbers are returned.
@@ -1006,17 +1008,9 @@ class CameraGeometry:
                 )
                 del dist_check
                 if index_check != insidepix_index:
-                    pix_indices[index] = -1
-
-        # print warning:
-        for index in np.where(pix_indices == -1)[0]:
-            logger.warning(
-                " Coordinate ({} m, {} m) lies outside camera".format(
-                    points_searched[0][index, 0], points_searched[0][index, 1]
-                )
-            )
+                    pix_indices[index] = invalid
 
-        return pix_indices if len(pix_indices) > 1 else pix_indices[0]
+        return np.squeeze(pix_indices) if scalar else pix_indices
 
     @staticmethod
     def simtel_shape_to_type(pixel_shape):
diff --git a/docs/changes/2397.api.rst b/docs/changes/2397.api.rst
new file mode 100644
index 00000000..45b64781
--- /dev/null
+++ b/docs/changes/2397.api.rst
@@ -0,0 +1,8 @@
+``CameraGeometry.position_to_pix_index`` will now return the minimum integer value for invalid
+pixel coordinates instead of -1 due to the danger of using -1 as an index in python accessing
+the last element of a data array for invalid pixels.
+The function will now also no longer raise an error if the arguments are empty arrays and instead
+just return an empty index array.
+The function will also no longer log a warning in case of coordinates that do not match a camera pixel.
+The function is very low-level and if not finding a pixel at the tested position warrants a warning or
+is expected will depend on the calling code.


Test patch:
diff --git a/ctapipe/instrument/camera/tests/test_geometry.py b/ctapipe/instrument/camera/tests/test_geometry.py
index 10410e0e..8e3a7545 100644
--- a/ctapipe/instrument/camera/tests/test_geometry.py
+++ b/ctapipe/instrument/camera/tests/test_geometry.py
@@ -68,10 +68,20 @@ def test_load_lst_camera(prod5_lst):
 
 def test_position_to_pix_index(prod5_lst):
     """test that we can lookup a pixel from a coordinate"""
+    geometry = prod5_lst.camera.geometry
+
     x, y = (0.80 * u.m, 0.79 * u.m)
-    pix_id = prod5_lst.camera.geometry.position_to_pix_index(x, y)
+
+    pix_id = geometry.position_to_pix_index(x, y)
+
     assert pix_id == 1575
 
+    pix_ids = geometry.position_to_pix_index([0.8, 0.8] * u.m, [0.79, 0.79] * u.m)
+    np.testing.assert_array_equal(pix_ids, [1575, 1575])
+
+    assert len(geometry.position_to_pix_index([] * u.m, [] * u.m)) == 0
+    assert geometry.position_to_pix_index(5 * u.m, 5 * u.m) == np.iinfo(int).min
+
 
 def test_find_neighbor_pixels():
     """test basic neighbor functionality"""