Repository: gammapy/gammapy (python)
Issue: FluxPoints.resample_axis fails in 1.3
**Gammapy version**
1.3, dev



**Bug description**
FluxPoints.resample_axis raises an error




**To Reproduce**
In the `analysis_2.py` tutorial, add at the end
```
from gammapy.estimators.utils import get_rebinned_axis
ax_new = get_rebinned_axis(fluxpoint=flux_points, method="min-ts", ts_threshold=10)
flux_points.resample_axis(ax_new)
```

This fails with an error:

```
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[40], line 1
----> 1 flux_points.resample_axis(ax_new)

File ~/Gammapy-dev/gammapy/gammapy/estimators/points/core.py:869, in FluxPoints.resample_axis(self, axis_new)
    867         edge_max = edge_max + axis_new.reference_time
    868     fp = self.slice_by_coord({axis_new.name: slice(edge_min, edge_max)})
--> 869     fp_new = squash_fluxpoints(fp, axis_new)
    870     fluxpoints.append(fp_new)
    872 return self.__class__.from_stack(fluxpoints, axis=axis_new)

File ~/Gammapy-dev/gammapy/gammapy/estimators/points/core.py:55, in squash_fluxpoints(flux_point, axis)
     52 if axis.name != "energy":
     53     geom = geom.to_cube([flux_point.geom.axes["energy"]])
---> 55 maps["norm"] = Map.from_geom(geom, data=minimizer.x)
     56 maps["norm_err"] = Map.from_geom(
     57     geom, data=np.sqrt(minimizer.hess_inv.todense()).reshape(geom.data_shape)
     58 )
     59 maps["n_dof"] = Map.from_geom(geom, data=flux_point.geom.axes[axis.name].nbin)

File ~/Gammapy-dev/gammapy/gammapy/maps/core.py:304, in Map.from_geom(geom, meta, data, unit, dtype)
    301     raise ValueError("Unrecognized geom type.")
    303 cls_out = Map._get_map_cls(map_type)
--> 304 return cls_out(geom, data=data, meta=meta, unit=unit, dtype=dtype)

File ~/Gammapy-dev/gammapy/gammapy/maps/region/ndmap.py:56, in RegionNDMap.__init__(self, geom, data, dtype, meta, unit)
     53     meta = {}
     55 self._geom = geom
---> 56 self.data = data
     57 self.meta = meta
     58 self._unit = u.Unit(unit)

File ~/Gammapy-dev/gammapy/gammapy/maps/region/ndmap.py:84, in RegionNDMap.data(self, value)
     81         value = np.expand_dims(value, (-2, -1))
     83 if self.geom.data_shape != value.shape:
---> 84     raise ValueError(
     85         f"Input shape {input_shape} is not compatible with shape from geometry {self.geom.data_shape}"
     86     )
     88 self._data = value

ValueError: Input shape (1,) is not compatible with shape from geometry (1, 1)
```



Pull request description: Fix bug in squash_fluxpoints
This fixes #5670 
The error was occurring because of the prevention of automatic reshaping introduced in #5180

Patch:
diff --git a/gammapy/estimators/points/core.py b/gammapy/estimators/points/core.py
index 58213992b..9407e5e97 100644
--- a/gammapy/estimators/points/core.py
+++ b/gammapy/estimators/points/core.py
@@ -52,32 +52,40 @@ def squash_fluxpoints(flux_point, axis):
     if axis.name != "energy":
         geom = geom.to_cube([flux_point.geom.axes["energy"]])
 
-    maps["norm"] = Map.from_geom(geom, data=minimizer.x)
+    maps["norm"] = Map.from_geom(geom, data=minimizer.x.reshape(geom.data_shape))
     maps["norm_err"] = Map.from_geom(
         geom, data=np.sqrt(minimizer.hess_inv.todense()).reshape(geom.data_shape)
     )
-    maps["n_dof"] = Map.from_geom(geom, data=flux_point.geom.axes[axis.name].nbin)
+    maps["n_dof"] = Map.from_geom(
+        geom, data=np.reshape(flux_point.geom.axes[axis.name].nbin, geom.data_shape)
+    )
 
     if "norm_ul" in flux_point.available_quantities:
         delta_ts = flux_point.meta.get("n_sigma_ul", 2) ** 2
         ul = stat_profile_ul_scipy(value_scan, stat_scan, delta_ts=delta_ts)
-        maps["norm_ul"] = Map.from_geom(geom, data=ul.value)
+        maps["norm_ul"] = Map.from_geom(
+            geom, data=np.reshape(ul.value, geom.data_shape)
+        )
 
-    maps["stat"] = Map.from_geom(geom, data=f(minimizer.x))
+    maps["stat"] = Map.from_geom(geom, data=f(minimizer.x).reshape(geom.data_shape))
 
     geom_scan = geom.to_cube([MapAxis.from_nodes(value_scan, name="norm")])
     maps["stat_scan"] = Map.from_geom(
         geom=geom_scan, data=stat_scan.reshape(geom_scan.data_shape)
     )
     try:
-        maps["stat_null"] = Map.from_geom(geom, data=np.sum(flux_point.stat_null.data))
+        maps["stat_null"] = Map.from_geom(
+            geom, data=np.reshape(np.sum(flux_point.stat_null.data), geom.data_shape)
+        )
         maps["ts"] = maps["stat_null"] - maps["stat"]
     except AttributeError:
         log.info(
             "Stat null info not present on original FluxPoints object. TS not computed"
         )
 
-    maps["success"] = Map.from_geom(geom=geom, data=minimizer.success, dtype=bool)
+    maps["success"] = Map.from_geom(
+        geom=geom, data=np.reshape(minimizer.success, geom.data_shape), dtype=bool
+    )
 
     combined_fp = FluxPoints.from_maps(
         maps=maps,


Test patch:
diff --git a/gammapy/estimators/points/tests/test_sed.py b/gammapy/estimators/points/tests/test_sed.py
index 1fb952efe..77f96e030 100644
--- a/gammapy/estimators/points/tests/test_sed.py
+++ b/gammapy/estimators/points/tests/test_sed.py
@@ -15,6 +15,7 @@ from gammapy.datasets import (
 )
 from gammapy.datasets.spectrum import SpectrumDataset
 from gammapy.estimators import FluxPoints, FluxPointsEstimator
+from gammapy.estimators.utils import get_rebinned_axis
 from gammapy.irf import EDispKernelMap, EffectiveAreaTable2D, load_irf_dict_from_file
 from gammapy.makers import MapDatasetMaker
 from gammapy.makers.utils import make_map_exposure_true_energy
@@ -853,3 +854,11 @@ def test_fpe_diff_lengths():
     datasets = Datasets([dataset1, dataset2, dataset3, dataset4])
     with pytest.raises(ValueError):
         fp = fpe.run(datasets)
+
+
+def test_resample_axis():
+    ds, fpe = create_fpe(PowerLawSpectralModel())
+    flux_points = fpe.run(ds)
+    axis_new = get_rebinned_axis(flux_points, method="fixed-bins", group_size=3)
+    flux_new = flux_points.resample_axis(axis_new)
+    assert_allclose(flux_new.flux.data, [[[3.1050191 * 1e-12]]], rtol=1e-5)