skymap inconsistency

[jballet]

Describe the bug
We use the fermipy utilities in fermipy.skymap to work with HEALPix maps.
We noted that the get_map_values call did not return data in the same format depending on whether the input is a HEALPix or a standard WCS map.
More precisely, with a HEALPix map the output of get_map_values has shape (ne,np), where ne is the number of energy bins and np the number of points in the sky.
With a WCS map the output has shape (np,ne) (transposed).
Obviously the same call cannot be used seamlessly with both kinds of inputs.
This is with fermipy version 1.2.0+2.gfa69.

With fermipy version 1.2.2+44.g580e it fails outright on
File "/dsm/manip/glast/anaconda3/envs/fermi-2.2.50/lib/python3.11/site-packages/fermipy/skymap.py", line 381, in get_map_values
vals[~m] = np.nan
~~~~^^^^
IndexError: boolean index did not match indexed array along dimension 0; dimension is 3 but corresponding boolean dimension is 4

To Reproduce
Here is a small test harness allowing to reproduce the bug. Just extract, go to the directory and call ./testskymap.py

Expected behavior
I would like the output to be (ne,np) (as in HEALPix) for WCS maps as well.

Environment (please complete the following information):
CentOS Linux release 7.8.2003 (Core)

[ndilalla]

Thanks @jballet for reporting this and providing the test. I will take a look and keep you updated.