Repository: gammapy/gammapy (python)

Issue:
WcsNDmap fill_by_coord creating strange wavy patterns in spectra
**Gammapy version**
Dev version

**Bug description** 
Creating a WcsNDmap with `fill_by_coord`, for example a spectrum, doesn't give the same results as plt.hist() as shown below :

![Capture d’écran 2024-07-04 à 16 57 19](https://github.com/gammapy/gammapy/assets/8652898/8e39e0f8-c1c2-4c05-901a-df0a18d73be3)

**Expected behavior**
Here I generated a fake event list with regularly spaced energies where one would expect 2 counts per energy bin.
We see wavy structure to assigning events sometimes to the left bin or to the right energy bin.
plt.hist() has correctly 2 counts in every bin.

Note that in the MRE if one changes 
evt_energy = np.arange(1000.,1020,0.5)*u.GeV 
to
evt_energy = np.arange(1000.0001,1020,0.5)*u.GeV 

this produces a flat spectrum as expected.

Issue could be related to [pix_tuple_to_idx](https://github.com/gammapy/gammapy/blob/2829b45f3a8b8c4885cc8a3d9984e919c171bab5/gammapy/maps/geom.py#L31) where `Pixel coordinates are rounded to the closest integer value.` 

So this seems to be a rounding issue

**To Reproduce**
Here is a minimal example : 

```
import numpy as np
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord
from astropy import units as u
from gammapy.maps import WcsGeom, WcsNDMap, MapAxis

E_edges=np.arange(1000,1021,1) *u.GeV #define some edges
evt_energy = np.arange(1000.,1020,0.5)*u.GeV #a fake photon event energy list 
print(f'E_edges= {E_edges}')
print(f'E_evt= {evt_energy}')

N_evt = evt_energy.shape[0]
radec = SkyCoord([center.ra]*N_evt, [center.dec]*N_evt,frame='fk5',unit="deg" ) #fake position list
energy_axis = MapAxis.from_edges(E_edges,name='energy')

map = WcsNDMap.create(skydir=center, binsz=0.1*u.deg, width=0.1*u.deg, frame='icrs', axes=[energy_axis], proj='CAR')
map.fill_by_coord({'skycoord': radec, 'energy': evt_energy})

spec = map.data.sum(axis=(1,2))
plt.plot(energy_axis.center, spec, marker='o',label='WcsNDMap fill_by_coord') 
h=plt.hist(evt_energy, bins=E_edges,histtype='step',label='plt.hist')
plt.legend()
plt.show()
# both histogram don't look alike why !!??
```



Pull request:
Remove use of np.rint in pix_tuple_to_idx
This fixes #5360 by using `np.floor(pixel+0.5)` instead of `np.rint(pixel)` 

The issue with `np.rint` (and `np.round`) is that `For values exactly halfway between rounded decimal values, NumPy
rounds to the nearest even value. Thus 1.5 and 2.5 round to 2.0, -0.5 and 0.5 round to 0.0, etc.`

The implementation here is similar to the implementation in `astropy.wcs.world_to_array_index_values` , https://docs.astropy.org/en/stable/_modules/astropy/wcs/wcsapi/low_level_api.html#BaseLowLevelWCS.world_to_array_index_values


This solves the issue in #5360 - the new plot given in the MRE  now matches 
![Unknown](https://github.com/gammapy/gammapy/assets/32677370/d322280f-b157-4c73-9a17-f417183e47bd)

I am labelling this as a bug fix, but we can decide if we want to backport it. This did not affect any of the test results, so probably this issue is evident only in artificially constructed examples. 

Patch:
diff --git a/gammapy/maps/geom.py b/gammapy/maps/geom.py
index 6a5d90efa..66165b7b2 100644
--- a/gammapy/maps/geom.py
+++ b/gammapy/maps/geom.py
@@ -47,7 +47,7 @@ def pix_tuple_to_idx(pix):
             idx += [p]
         else:
             with np.errstate(invalid="ignore"):
-                p_idx = np.rint(p).astype(int)
+                p_idx = np.floor(p + 0.5).astype(int)
             p_idx[~np.isfinite(p)] = INVALID_INDEX.int
             idx += [p_idx]
 


Test patch:
diff --git a/gammapy/maps/wcs/tests/test_ndmap.py b/gammapy/maps/wcs/tests/test_ndmap.py
index 022ba4a99..9ed7a09b5 100644
--- a/gammapy/maps/wcs/tests/test_ndmap.py
+++ b/gammapy/maps/wcs/tests/test_ndmap.py
@@ -1045,3 +1045,26 @@ def test_plot_mask():
 
     with mpl_plot_check():
         mask.plot_grid()
+
+
+def test_histogram_center_value():
+    evt_energy = np.arange(100.0, 105, 0.5) * u.GeV
+    N_evt = evt_energy.shape[0]
+    center = SkyCoord(0, 0, unit="deg", frame="icrs")
+    radec = SkyCoord(
+        [center.ra] * N_evt, [center.dec] * N_evt, unit="deg"
+    )  # fake position list
+    energy_axis = MapAxis.from_edges(np.arange(100, 106, 1) * u.GeV, name="energy")
+    map1 = WcsNDMap.create(
+        skydir=center, binsz=0.1 * u.deg, width=0.1 * u.deg, axes=[energy_axis]
+    )
+    map1.fill_by_coord({"skycoord": radec, "energy": evt_energy})
+    assert_allclose(map1.data[0:2], [[[2.0]], [[2.0]]])
+
+    map1 = WcsNDMap.create(
+        skydir=center,
+        binsz=1 * u.deg,
+        width=3 * u.deg,
+    )
+    map1.fill_by_pix([[0.0, 0.5, 1.0, 1.5], [0.0, 0.5, 1.0, 1.5]])
+    assert_allclose(map1.data, [[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 1.0]])