WCSNDMap reorders input data

[maxnoe]

Gammapy version

Current main branch and also released versions.

Bug description

The WCSNDMap silently reorders the array data by calling data.reshape(geom.data_shape) if the input data shape does not match the expected shape.

This is most certainly wrong, as you cannot just reshape the same memory into the desired shape and assume the values are correct.

E.g. in the following example, the correct thing would probably have been a transpose operation (I assume WCSGeom assumes FITS / Fortran shape order), not the reshaping operation.

In [1]: from gammapy.maps import WcsNDMap, WcsGeom
   ...: import numpy as np
   ...: from astropy.coordinates import SkyCoord

In [2]: data = np.arange(15).reshape(5, 3)

In [3]: center = SkyCoord.from_name("Crab Nebula")

In [4]: geom = WcsGeom.create(
   ...:     npix=(5, 3),
   ...:     binsz=0.5,
   ...:     refpix=(2, 1),
   ...:     skydir=center,
   ...: )

In [5]: wcsmap = WcsNDMap(geom, data)

In [6]: wcsmap.data
Out[6]: 
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14]])

In [7]: data
Out[7]: 
array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11],
       [12, 13, 14]])

In [8]: wcsmap.data.shape == data.T.shape
Out[8]: True

In [9]: wcsmap.data == data.T
Out[9]: 
array([[ True, False, False, False, False],
       [False, False,  True, False, False],
       [False, False, False, False,  True]])

Expected behavior

An error in case the data shape does not match instead of silently corrupting the input data.

To Reproduce

See example above.

Other information

CC @StFroese @LukasNickel