Repository: gammapy/gammapy (python)

Issue:
WcsNDMap.downsample() `weights` parameter normalization issue
**Gammapy version**
1.2


**Bug description**
When downsampling a WcsNDMap with `preserve_counts=False`, i.e. averaging, the weighted average is not properly calculated. Instead an average of `weights.data*self.data` is returned.


**Expected behavior**
The expected formula for a weighted average would be `np.sum(weights.data*self.data)/np.sum(weights.data)`.


**To Reproduce**

Loading an example dataset from the tutorials:
```python
dataset_cta = MapDataset.read(
    "$GAMMAPY_DATA/cta-1dc-gc/cta-1dc-gc.fits.gz", name="dataset-cta")
```
Compare the scales of
```python
dataset_cta.psf.psf_map.plot_interactive(add_cbar=True)
```
and

```python
dataset_cta.psf.psf_map.downsample(
    2,preserve_counts=False,weights=dataset_cta.psf.exposure_map).plot_interactive(add_cbar=True)
```

The result should be a downsampled map weighted with the exposure map, the resulting scales are, however, around 1e9 too high, as the exposure map values have not been normalized out.


**Other information**
The downsampled map data is calculated in this line in the WcsNDMap definition:
```python
data = block_reduce(self.data * weights, tuple(block_size), func=func)
```
which works only for `preserve_counts=True`, where `func=np.nansum`, but not for `preserve_counts=False`, where `func=np.nanmean`.

A solution could be, for instance:
```python
data = block_reduce(self.data * weights, tuple(block_size), func=np.nansum)
if not preserve_counts:
    weight_sum = block_reduce(weights, tuple(block_size), func=np.nansum)
    data=np.divide(data,weight_sum, where=np.where(data>0.0,True,False))
```

Please tell me if I've missed anything here! Thank you!

Pull request description:
WcsNDMap.downsample normalization issue fix
<!-- These are hidden commments. You can leave them or delete them. -->

<!-- Don't be shy! Github issues or pull requests are welcome any time. -->
<!-- But if you first want to chat, use the Gammapy Slack. -->
<!-- See https://gammapy.org/contact.html -->

**Description**
<!-- What is the motivation for this pull request? -->
<!-- Briefly: what changes are done here? -->
<!-- If this is to address a Github issue, mention its number and use the specific keywords to create a link -->
<!-- See docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue -->
<!-- Example: This PR is to resolve #42 -->

This pull request fixes #5493. Instead of specifying either `np.nansum` or `np.nanmean` as functions for `block_reduce`, `self.data * weights` is now always summed and afterwards divided by the sum of the weights. This way the weighted average is calculated correctly.

**Dear reviewer**
<!-- Let the reviewer and Gammapy team know what you want: -->
<!-- * Is this ready for review? Or is it work in progress and you want some feedback? -->
<!-- * Do you want to go through review here? Or if someone should just finish this up and merge it? -->
<!-- Do you have any specific questions, e.g. about API or implementation? -->
<!-- Do you include a test executing the new code you're adding (to make sure it runs)? -->
<!-- Do you include some documentation? Is it needed? -->

This pull request is ready for review! Please let me know if anything should be changed or I missed anything. Thank you!


Patch:
diff --git a/gammapy/maps/wcs/ndmap.py b/gammapy/maps/wcs/ndmap.py
index 97f34e237..02946beea 100644
--- a/gammapy/maps/wcs/ndmap.py
+++ b/gammapy/maps/wcs/ndmap.py
@@ -358,14 +358,22 @@ class WcsNDMap(WcsMap):
             idx = self.geom.axes.index_data(axis_name)
             block_size[idx] = factor
 
-        func = np.nansum if preserve_counts else np.nanmean
-
         if weights is None:
-            weights = 1
-        else:
+            weights = np.ones_like(self.data)
+        elif not isinstance(weights, np.ndarray):
             weights = weights.data
 
-        data = block_reduce(self.data * weights, tuple(block_size), func=func)
+        data = block_reduce(self.data * weights, tuple(block_size), func=np.nansum)
+
+        if not preserve_counts:
+            weight_sum = block_reduce(weights, tuple(block_size), func=np.nansum)
+            data = np.divide(
+                data,
+                weight_sum,
+                out=np.zeros_like(data, dtype="float64"),
+                where=(weight_sum != 0),
+            )
+
         return self._init_copy(geom=geom, data=data.astype(self.data.dtype))
 
     def plot(


Test patch:
diff --git a/gammapy/maps/wcs/tests/test_ndmap.py b/gammapy/maps/wcs/tests/test_ndmap.py
index 1101edb60..bd34cc301 100644
--- a/gammapy/maps/wcs/tests/test_ndmap.py
+++ b/gammapy/maps/wcs/tests/test_ndmap.py
@@ -426,10 +426,42 @@ def test_wcsndmap_downsample(npix, binsz, frame, proj, skydir, axes):
     m = WcsNDMap(geom, unit="m2")
     # Check whether we can downsample
     if np.all(np.mod(geom.npix[0], 2) == 0) and np.all(np.mod(geom.npix[1], 2) == 0):
+        # without weights
         m2 = m.downsample(2, preserve_counts=True)
         assert_allclose(np.nansum(m.data), np.nansum(m2.data))
         assert m.unit == m2.unit
 
+        m3 = m.downsample(2, preserve_counts=False)
+        assert_allclose(np.nanmean(m.data), np.nanmean(m3.data))
+        assert m.unit == m3.unit
+
+        # with weights
+        weights = np.random.rand(*m.data.shape)
+        m4 = m.downsample(2, preserve_counts=True, weights=weights)
+        assert_allclose(np.nansum(m.data * weights.data), np.nansum(m4.data))
+        assert m.unit == m4.unit
+
+
+def test_wcsndmap_downsample_weights():
+    geom = WcsGeom.create(
+        npix=10,
+        binsz=1.0,
+        frame="galactic",
+        proj="AIT",
+        skydir=SkyCoord(110.0, 75.0, unit="deg", frame="icrs"),
+        axes=None,
+    )
+    m = WcsNDMap(geom, data=np.tile(np.array([3.0, 2.0]), (10, 5)), unit="m2")
+
+    # weights array with recurring values to easily test for correct weighted averaging
+    weights = np.tile(np.array([1.0, 2.0]), (m.data.shape[0], m.data.shape[1] // 2))
+    m6 = m.downsample(2, preserve_counts=False, weights=weights)
+
+    assert_allclose(
+        np.nansum(m.data * weights) / np.nansum(weights), np.nanmean(m6.data)
+    )
+    assert m.unit == m6.unit
+
 
 @pytest.mark.parametrize(
     ("npix", "binsz", "frame", "proj", "skydir", "axes"), wcs_test_geoms