We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Fix test_wcs_nd_map_data_transpose_issue
I think here we forgot to pass `sparse=True` in `test_wcs_nd_map_data_transpose_issue`?

https://github.com/gammapy/gammapy/pull/1346/files#diff-930fa0f94b49b5d6537396fad7c6512aR137-R140

If I change to `sparse=True` in the write call, I get this fail:
```
E           Arrays are not equal
E           
E           x and y nan location mismatch:
E            x: array([[0., 1., 2.],
E                  [0., 0., 0.]], dtype=float32)
E            y: array([[  0.,   1.,   2.],
E                  [ nan,  inf, -inf]])
```

@adonath - Is there a bug?

Should we spend 1-2 days on maps basics and sort out dtype and missing value issues (#1440, #2374, maybe this issue), e.g. with pair-coding & debugging? Or defer all or some to v1.1?

This is a reminder issue, noticed this accidentally while working on
https://github.com/gammapy/gammapy/pull/2350#issuecomment-542754564

--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `gammapy/maps/wcsmap.py`
Content:
```
1 # Licensed under a 3-clause BSD style license - see LICENSE.rst
2 import json
3 import numpy as np
4 from astropy.io import fits
5 from .core import Map
6 from .utils import find_bands_hdu, find_hdu, JsonQuantityEncoder
7 from .wcs import WcsGeom
8 
9 
10 __all__ = ["WcsMap"]
11 
12 
13 def identify_wcs_format(hdu):
14     if hdu is None:
15         return "gadf"
16     elif hdu.name == "ENERGIES":
17         return "fgst-template"
18     elif hdu.name == "EBOUNDS":
19         return "fgst-ccube"
20     else:
21         return "gadf"
22 
23 
24 class WcsMap(Map):
25     """Base class for WCS map classes.
26 
27     Parameters
28     ----------
29     geom : `~gammapy.maps.WcsGeom`
30         A WCS geometry object.
31     data : `~numpy.ndarray`
32         Data array.
33     """
34 
35     @classmethod
36     def create(
37         cls,
38         map_type="wcs",
39         npix=None,
40         binsz=0.1,
41         width=None,
42         proj="CAR",
43         frame="icrs",
44         refpix=None,
45         axes=None,
46         skydir=None,
47         dtype="float32",
48         meta=None,
49         unit="",
50     ):
51         """Factory method to create an empty WCS map.
52 
53         Parameters
54         ----------
55         map_type : {'wcs', 'wcs-sparse'}
56             Map type.  Selects the class that will be used to
57             instantiate the map.
58         npix : int or tuple or list
59             Width of the map in pixels. A tuple will be interpreted as
60             parameters for longitude and latitude axes.  For maps with
61             non-spatial dimensions, list input can be used to define a
62             different map width in each image plane.  This option
63             supersedes width.
64         width : float or tuple or list
65             Width of the map in degrees.  A tuple will be interpreted
66             as parameters for longitude and latitude axes.  For maps
67             with non-spatial dimensions, list input can be used to
68             define a different map width in each image plane.
69         binsz : float or tuple or list
70             Map pixel size in degrees.  A tuple will be interpreted
71             as parameters for longitude and latitude axes.  For maps
72             with non-spatial dimensions, list input can be used to
73             define a different bin size in each image plane.
74         skydir : tuple or `~astropy.coordinates.SkyCoord`
75             Sky position of map center.  Can be either a SkyCoord
76             object or a tuple of longitude and latitude in deg in the
77             coordinate system of the map.
78         frame : {"icrs", "galactic"}, optional
79             Coordinate system, either Galactic ("galactic") or Equatorial ("icrs").
80         axes : list
81             List of non-spatial axes.
82         proj : string, optional
83             Any valid WCS projection type. Default is 'CAR' (cartesian).
84         refpix : tuple
85             Reference pixel of the projection.  If None then this will
86             be chosen to be center of the map.
87         dtype : str, optional
88             Data type, default is float32
89         meta : `dict`
90             Dictionary to store meta data.
91         unit : str or `~astropy.units.Unit`
92             The unit of the map
93 
94         Returns
95         -------
96         map : `~WcsMap`
97             A WCS map object.
98         """
99         from .wcsnd import WcsNDMap
100 
101         geom = WcsGeom.create(
102             npix=npix,
103             binsz=binsz,
104             width=width,
105             proj=proj,
106             skydir=skydir,
107             frame=frame,
108             refpix=refpix,
109             axes=axes,
110         )
111 
112         if map_type == "wcs":
113             return WcsNDMap(geom, dtype=dtype, meta=meta, unit=unit)
114         elif map_type == "wcs-sparse":
115             raise NotImplementedError
116         else:
117             raise ValueError(f"Invalid map type: {map_type!r}")
118 
119     @classmethod
120     def from_hdulist(cls, hdu_list, hdu=None, hdu_bands=None, format="gadf"):
121         """Make a WcsMap object from a FITS HDUList.
122 
123         Parameters
124         ----------
125         hdu_list :  `~astropy.io.fits.HDUList`
126             HDU list containing HDUs for map data and bands.
127         hdu : str
128             Name or index of the HDU with the map data.
129         hdu_bands : str
130             Name or index of the HDU with the BANDS table.
131         format : {'gadf', 'fgst-ccube', 'fgst-template'}
132             FITS format convention.
133 
134         Returns
135         -------
136         wcs_map : `WcsMap`
137             Map object
138         """
139         if hdu is None:
140             hdu = find_hdu(hdu_list)
141         else:
142             hdu = hdu_list[hdu]
143 
144         if hdu_bands is None:
145             hdu_bands = find_bands_hdu(hdu_list, hdu)
146 
147         if hdu_bands is not None:
148             hdu_bands = hdu_list[hdu_bands]
149 
150         format = identify_wcs_format(hdu_bands)
151 
152         wcs_map = cls.from_hdu(hdu, hdu_bands, format=format)
153 
154         if wcs_map.unit.is_equivalent(""):
155             if format == "fgst-template":
156                 if "GTI" in hdu_list:  # exposure maps have an additional GTI hdu
157                     wcs_map.unit = "cm2 s"
158                 else:
159                     wcs_map.unit = "cm-2 s-1 MeV-1 sr-1"
160 
161         return wcs_map
162 
163     def to_hdulist(self, hdu=None, hdu_bands=None, sparse=False, format="gadf"):
164         """Convert to `~astropy.io.fits.HDUList`.
165 
166         Parameters
167         ----------
168         hdu : str
169             Name or index of the HDU with the map data.
170         hdu_bands : str
171             Name or index of the HDU with the BANDS table.
172         sparse : bool
173             Sparsify the map by only writing pixels with non-zero
174             amplitude.
175         format : {'gadf', 'fgst-ccube','fgst-template'}
176             FITS format convention.
177 
178         Returns
179         -------
180         hdu_list : `~astropy.io.fits.HDUList`
181 
182         """
183         if sparse:
184             hdu = "SKYMAP" if hdu is None else hdu.upper()
185         else:
186             hdu = "PRIMARY" if hdu is None else hdu.upper()
187 
188         if sparse and hdu == "PRIMARY":
189             raise ValueError("Sparse maps cannot be written to the PRIMARY HDU.")
190 
191         if format in ["fgst-ccube", "fgst-template"]:
192             if self.geom.axes[0].name != "energy" or len(self.geom.axes) > 1:
193                 raise ValueError(
194                     "All 'fgst' formats don't support extra axes except for energy."
195                 )
196 
197         if hdu_bands is None:
198             hdu_bands = f"{hdu.upper()}_BANDS"
199 
200         if self.geom.axes:
201             hdu_bands_out = self.geom.to_bands_hdu(
202                 hdu_bands=hdu_bands, format=format
203             )
204             hdu_bands = hdu_bands_out.name
205         else:
206             hdu_bands = None
207 
208         hdu_out = self.to_hdu(hdu=hdu, hdu_bands=hdu_bands, sparse=sparse)
209 
210         hdu_out.header["META"] = json.dumps(self.meta, cls=JsonQuantityEncoder)
211 
212         hdu_out.header["BUNIT"] = self.unit.to_string("fits")
213 
214         if hdu == "PRIMARY":
215             hdulist = [hdu_out]
216         else:
217             hdulist = [fits.PrimaryHDU(), hdu_out]
218 
219         if self.geom.axes:
220             hdulist += [hdu_bands_out]
221 
222         return fits.HDUList(hdulist)
223 
224     def to_hdu(self, hdu="SKYMAP", hdu_bands=None, sparse=False):
225         """Make a FITS HDU from this map.
226 
227         Parameters
228         ----------
229         hdu : str
230             The HDU extension name.
231         hdu_bands : str
232             The HDU extension name for BANDS table.
233         sparse : bool
234             Set INDXSCHM to SPARSE and sparsify the map by only
235             writing pixels with non-zero amplitude.
236 
237         Returns
238         -------
239         hdu : `~astropy.io.fits.BinTableHDU` or `~astropy.io.fits.ImageHDU`
240             HDU containing the map data.
241         """
242         header = self.geom.to_header()
243 
244         if self.is_mask:
245             data = self.data.astype(int)
246         else:
247             data = self.data
248 
249         if hdu_bands is not None:
250             header["BANDSHDU"] = hdu_bands
251 
252         if sparse:
253             hdu_out = self._make_hdu_sparse(data, self.geom.npix, hdu, header)
254         elif hdu == "PRIMARY":
255             hdu_out = fits.PrimaryHDU(data, header=header)
256         else:
257             hdu_out = fits.ImageHDU(data, header=header, name=hdu)
258 
259         return hdu_out
260 
261     @staticmethod
262     def _make_hdu_sparse(data, npix, hdu, header):
263         shape = data.shape
264 
265         # We make a copy, because below we modify `data` to handle non-finite entries
266         # TODO: The code below could probably be simplified to use expressions
267         # that create new arrays instead of in-place modifications
268         # But first: do we want / need the non-finite entry handling at all and always cast to 64-bit float?
269         data = data.copy()
270 
271         if len(shape) == 2:
272             data_flat = np.ravel(data)
273             data_flat[~np.isfinite(data_flat)] = 0
274             nonzero = np.where(data_flat > 0)
275             value = data_flat[nonzero].astype(float)
276             cols = [
277                 fits.Column("PIX", "J", array=nonzero[0]),
278                 fits.Column("VALUE", "E", array=value),
279             ]
280         elif npix[0].size == 1:
281             shape_flat = shape[:-2] + (shape[-1] * shape[-2],)
282             data_flat = np.ravel(data).reshape(shape_flat)
283             data_flat[~np.isfinite(data_flat)] = 0
284             nonzero = np.where(data_flat > 0)
285             channel = np.ravel_multi_index(nonzero[:-1], shape[:-2])
286             value = data_flat[nonzero].astype(float)
287             cols = [
288                 fits.Column("PIX", "J", array=nonzero[-1]),
289                 fits.Column("CHANNEL", "I", array=channel),
290                 fits.Column("VALUE", "E", array=value),
291             ]
292         else:
293             data_flat = []
294             channel = []
295             pix = []
296             for i, _ in np.ndenumerate(npix[0]):
297                 data_i = np.ravel(data[i[::-1]])
298                 data_i[~np.isfinite(data_i)] = 0
299                 pix_i = np.where(data_i > 0)
300                 data_i = data_i[pix_i]
301                 data_flat += [data_i]
302                 pix += pix_i
303                 channel += [
304                     np.ones(data_i.size, dtype=int)
305                     * np.ravel_multi_index(i[::-1], shape[:-2])
306                 ]
307 
308             pix = np.concatenate(pix)
309             channel = np.concatenate(channel)
310             value = np.concatenate(data_flat).astype(float)
311 
312             cols = [
313                 fits.Column("PIX", "J", array=pix),
314                 fits.Column("CHANNEL", "I", array=channel),
315                 fits.Column("VALUE", "E", array=value),
316             ]
317 
318         return fits.BinTableHDU.from_columns(cols, header=header, name=hdu)
319
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/gammapy/maps/wcsmap.py b/gammapy/maps/wcsmap.py
--- a/gammapy/maps/wcsmap.py
+++ b/gammapy/maps/wcsmap.py
@@ -270,18 +270,16 @@
 
         if len(shape) == 2:
             data_flat = np.ravel(data)
-            data_flat[~np.isfinite(data_flat)] = 0
-            nonzero = np.where(data_flat > 0)
-            value = data_flat[nonzero].astype(float)
+            non_zero = np.where(~(data_flat == 0))
+            value = data_flat[non_zero].astype(float)
             cols = [
-                fits.Column("PIX", "J", array=nonzero[0]),
+                fits.Column("PIX", "J", array=non_zero[0]),
                 fits.Column("VALUE", "E", array=value),
             ]
         elif npix[0].size == 1:
             shape_flat = shape[:-2] + (shape[-1] * shape[-2],)
             data_flat = np.ravel(data).reshape(shape_flat)
-            data_flat[~np.isfinite(data_flat)] = 0
-            nonzero = np.where(data_flat > 0)
+            nonzero = np.where(~(data_flat == 0))
             channel = np.ravel_multi_index(nonzero[:-1], shape[:-2])
             value = data_flat[nonzero].astype(float)
             cols = [
@@ -295,8 +293,7 @@
             pix = []
             for i, _ in np.ndenumerate(npix[0]):
                 data_i = np.ravel(data[i[::-1]])
-                data_i[~np.isfinite(data_i)] = 0
-                pix_i = np.where(data_i > 0)
+                pix_i = np.where(~(data_i == 0))
                 data_i = data_i[pix_i]
                 data_flat += [data_i]
                 pix += pix_i

{"golden_diff": "diff --git a/gammapy/maps/wcsmap.py b/gammapy/maps/wcsmap.py\n--- a/gammapy/maps/wcsmap.py\n+++ b/gammapy/maps/wcsmap.py\n@@ -270,18 +270,16 @@\n \n         if len(shape) == 2:\n             data_flat = np.ravel(data)\n-            data_flat[~np.isfinite(data_flat)] = 0\n-            nonzero = np.where(data_flat > 0)\n-            value = data_flat[nonzero].astype(float)\n+            non_zero = np.where(~(data_flat == 0))\n+            value = data_flat[non_zero].astype(float)\n             cols = [\n-                fits.Column(\"PIX\", \"J\", array=nonzero[0]),\n+                fits.Column(\"PIX\", \"J\", array=non_zero[0]),\n                 fits.Column(\"VALUE\", \"E\", array=value),\n             ]\n         elif npix[0].size == 1:\n             shape_flat = shape[:-2] + (shape[-1] * shape[-2],)\n             data_flat = np.ravel(data).reshape(shape_flat)\n-            data_flat[~np.isfinite(data_flat)] = 0\n-            nonzero = np.where(data_flat > 0)\n+            nonzero = np.where(~(data_flat == 0))\n             channel = np.ravel_multi_index(nonzero[:-1], shape[:-2])\n             value = data_flat[nonzero].astype(float)\n             cols = [\n@@ -295,8 +293,7 @@\n             pix = []\n             for i, _ in np.ndenumerate(npix[0]):\n                 data_i = np.ravel(data[i[::-1]])\n-                data_i[~np.isfinite(data_i)] = 0\n-                pix_i = np.where(data_i > 0)\n+                pix_i = np.where(~(data_i == 0))\n                 data_i = data_i[pix_i]\n                 data_flat += [data_i]\n                 pix += pix_i\n", "issue": "Fix test_wcs_nd_map_data_transpose_issue\nI think here we forgot to pass `sparse=True` in `test_wcs_nd_map_data_transpose_issue`?\r\n\r\nhttps://github.com/gammapy/gammapy/pull/1346/files#diff-930fa0f94b49b5d6537396fad7c6512aR137-R140\r\n\r\nIf I change to `sparse=True` in the write call, I get this fail:\r\n```\r\nE           Arrays are not equal\r\nE           \r\nE           x and y nan location mismatch:\r\nE            x: array([[0., 1., 2.],\r\nE                  [0., 0., 0.]], dtype=float32)\r\nE            y: array([[  0.,   1.,   2.],\r\nE                  [ nan,  inf, -inf]])\r\n```\r\n\r\n@adonath - Is there a bug?\r\n\r\nShould we spend 1-2 days on maps basics and sort out dtype and missing value issues (#1440, #2374, maybe this issue), e.g. with pair-coding & debugging? Or defer all or some to v1.1?\r\n\r\nThis is a reminder issue, noticed this accidentally while working on\r\nhttps://github.com/gammapy/gammapy/pull/2350#issuecomment-542754564\n", "before_files": [{"content": "# Licensed under a 3-clause BSD style license - see LICENSE.rst\nimport json\nimport numpy as np\nfrom astropy.io import fits\nfrom .core import Map\nfrom .utils import find_bands_hdu, find_hdu, JsonQuantityEncoder\nfrom .wcs import WcsGeom\n\n\n__all__ = [\"WcsMap\"]\n\n\ndef identify_wcs_format(hdu):\n    if hdu is None:\n        return \"gadf\"\n    elif hdu.name == \"ENERGIES\":\n        return \"fgst-template\"\n    elif hdu.name == \"EBOUNDS\":\n        return \"fgst-ccube\"\n    else:\n        return \"gadf\"\n\n\nclass WcsMap(Map):\n    \"\"\"Base class for WCS map classes.\n\n    Parameters\n    ----------\n    geom : `~gammapy.maps.WcsGeom`\n        A WCS geometry object.\n    data : `~numpy.ndarray`\n        Data array.\n    \"\"\"\n\n    @classmethod\n    def create(\n        cls,\n        map_type=\"wcs\",\n        npix=None,\n        binsz=0.1,\n        width=None,\n        proj=\"CAR\",\n        frame=\"icrs\",\n        refpix=None,\n        axes=None,\n        skydir=None,\n        dtype=\"float32\",\n        meta=None,\n        unit=\"\",\n    ):\n        \"\"\"Factory method to create an empty WCS map.\n\n        Parameters\n        ----------\n        map_type : {'wcs', 'wcs-sparse'}\n            Map type.  Selects the class that will be used to\n            instantiate the map.\n        npix : int or tuple or list\n            Width of the map in pixels. A tuple will be interpreted as\n            parameters for longitude and latitude axes.  For maps with\n            non-spatial dimensions, list input can be used to define a\n            different map width in each image plane.  This option\n            supersedes width.\n        width : float or tuple or list\n            Width of the map in degrees.  A tuple will be interpreted\n            as parameters for longitude and latitude axes.  For maps\n            with non-spatial dimensions, list input can be used to\n            define a different map width in each image plane.\n        binsz : float or tuple or list\n            Map pixel size in degrees.  A tuple will be interpreted\n            as parameters for longitude and latitude axes.  For maps\n            with non-spatial dimensions, list input can be used to\n            define a different bin size in each image plane.\n        skydir : tuple or `~astropy.coordinates.SkyCoord`\n            Sky position of map center.  Can be either a SkyCoord\n            object or a tuple of longitude and latitude in deg in the\n            coordinate system of the map.\n        frame : {\"icrs\", \"galactic\"}, optional\n            Coordinate system, either Galactic (\"galactic\") or Equatorial (\"icrs\").\n        axes : list\n            List of non-spatial axes.\n        proj : string, optional\n            Any valid WCS projection type. Default is 'CAR' (cartesian).\n        refpix : tuple\n            Reference pixel of the projection.  If None then this will\n            be chosen to be center of the map.\n        dtype : str, optional\n            Data type, default is float32\n        meta : `dict`\n            Dictionary to store meta data.\n        unit : str or `~astropy.units.Unit`\n            The unit of the map\n\n        Returns\n        -------\n        map : `~WcsMap`\n            A WCS map object.\n        \"\"\"\n        from .wcsnd import WcsNDMap\n\n        geom = WcsGeom.create(\n            npix=npix,\n            binsz=binsz,\n            width=width,\n            proj=proj,\n            skydir=skydir,\n            frame=frame,\n            refpix=refpix,\n            axes=axes,\n        )\n\n        if map_type == \"wcs\":\n            return WcsNDMap(geom, dtype=dtype, meta=meta, unit=unit)\n        elif map_type == \"wcs-sparse\":\n            raise NotImplementedError\n        else:\n            raise ValueError(f\"Invalid map type: {map_type!r}\")\n\n    @classmethod\n    def from_hdulist(cls, hdu_list, hdu=None, hdu_bands=None, format=\"gadf\"):\n        \"\"\"Make a WcsMap object from a FITS HDUList.\n\n        Parameters\n        ----------\n        hdu_list :  `~astropy.io.fits.HDUList`\n            HDU list containing HDUs for map data and bands.\n        hdu : str\n            Name or index of the HDU with the map data.\n        hdu_bands : str\n            Name or index of the HDU with the BANDS table.\n        format : {'gadf', 'fgst-ccube', 'fgst-template'}\n            FITS format convention.\n\n        Returns\n        -------\n        wcs_map : `WcsMap`\n            Map object\n        \"\"\"\n        if hdu is None:\n            hdu = find_hdu(hdu_list)\n        else:\n            hdu = hdu_list[hdu]\n\n        if hdu_bands is None:\n            hdu_bands = find_bands_hdu(hdu_list, hdu)\n\n        if hdu_bands is not None:\n            hdu_bands = hdu_list[hdu_bands]\n\n        format = identify_wcs_format(hdu_bands)\n\n        wcs_map = cls.from_hdu(hdu, hdu_bands, format=format)\n\n        if wcs_map.unit.is_equivalent(\"\"):\n            if format == \"fgst-template\":\n                if \"GTI\" in hdu_list:  # exposure maps have an additional GTI hdu\n                    wcs_map.unit = \"cm2 s\"\n                else:\n                    wcs_map.unit = \"cm-2 s-1 MeV-1 sr-1\"\n\n        return wcs_map\n\n    def to_hdulist(self, hdu=None, hdu_bands=None, sparse=False, format=\"gadf\"):\n        \"\"\"Convert to `~astropy.io.fits.HDUList`.\n\n        Parameters\n        ----------\n        hdu : str\n            Name or index of the HDU with the map data.\n        hdu_bands : str\n            Name or index of the HDU with the BANDS table.\n        sparse : bool\n            Sparsify the map by only writing pixels with non-zero\n            amplitude.\n        format : {'gadf', 'fgst-ccube','fgst-template'}\n            FITS format convention.\n\n        Returns\n        -------\n        hdu_list : `~astropy.io.fits.HDUList`\n\n        \"\"\"\n        if sparse:\n            hdu = \"SKYMAP\" if hdu is None else hdu.upper()\n        else:\n            hdu = \"PRIMARY\" if hdu is None else hdu.upper()\n\n        if sparse and hdu == \"PRIMARY\":\n            raise ValueError(\"Sparse maps cannot be written to the PRIMARY HDU.\")\n\n        if format in [\"fgst-ccube\", \"fgst-template\"]:\n            if self.geom.axes[0].name != \"energy\" or len(self.geom.axes) > 1:\n                raise ValueError(\n                    \"All 'fgst' formats don't support extra axes except for energy.\"\n                )\n\n        if hdu_bands is None:\n            hdu_bands = f\"{hdu.upper()}_BANDS\"\n\n        if self.geom.axes:\n            hdu_bands_out = self.geom.to_bands_hdu(\n                hdu_bands=hdu_bands, format=format\n            )\n            hdu_bands = hdu_bands_out.name\n        else:\n            hdu_bands = None\n\n        hdu_out = self.to_hdu(hdu=hdu, hdu_bands=hdu_bands, sparse=sparse)\n\n        hdu_out.header[\"META\"] = json.dumps(self.meta, cls=JsonQuantityEncoder)\n\n        hdu_out.header[\"BUNIT\"] = self.unit.to_string(\"fits\")\n\n        if hdu == \"PRIMARY\":\n            hdulist = [hdu_out]\n        else:\n            hdulist = [fits.PrimaryHDU(), hdu_out]\n\n        if self.geom.axes:\n            hdulist += [hdu_bands_out]\n\n        return fits.HDUList(hdulist)\n\n    def to_hdu(self, hdu=\"SKYMAP\", hdu_bands=None, sparse=False):\n        \"\"\"Make a FITS HDU from this map.\n\n        Parameters\n        ----------\n        hdu : str\n            The HDU extension name.\n        hdu_bands : str\n            The HDU extension name for BANDS table.\n        sparse : bool\n            Set INDXSCHM to SPARSE and sparsify the map by only\n            writing pixels with non-zero amplitude.\n\n        Returns\n        -------\n        hdu : `~astropy.io.fits.BinTableHDU` or `~astropy.io.fits.ImageHDU`\n            HDU containing the map data.\n        \"\"\"\n        header = self.geom.to_header()\n\n        if self.is_mask:\n            data = self.data.astype(int)\n        else:\n            data = self.data\n\n        if hdu_bands is not None:\n            header[\"BANDSHDU\"] = hdu_bands\n\n        if sparse:\n            hdu_out = self._make_hdu_sparse(data, self.geom.npix, hdu, header)\n        elif hdu == \"PRIMARY\":\n            hdu_out = fits.PrimaryHDU(data, header=header)\n        else:\n            hdu_out = fits.ImageHDU(data, header=header, name=hdu)\n\n        return hdu_out\n\n    @staticmethod\n    def _make_hdu_sparse(data, npix, hdu, header):\n        shape = data.shape\n\n        # We make a copy, because below we modify `data` to handle non-finite entries\n        # TODO: The code below could probably be simplified to use expressions\n        # that create new arrays instead of in-place modifications\n        # But first: do we want / need the non-finite entry handling at all and always cast to 64-bit float?\n        data = data.copy()\n\n        if len(shape) == 2:\n            data_flat = np.ravel(data)\n            data_flat[~np.isfinite(data_flat)] = 0\n            nonzero = np.where(data_flat > 0)\n            value = data_flat[nonzero].astype(float)\n            cols = [\n                fits.Column(\"PIX\", \"J\", array=nonzero[0]),\n                fits.Column(\"VALUE\", \"E\", array=value),\n            ]\n        elif npix[0].size == 1:\n            shape_flat = shape[:-2] + (shape[-1] * shape[-2],)\n            data_flat = np.ravel(data).reshape(shape_flat)\n            data_flat[~np.isfinite(data_flat)] = 0\n            nonzero = np.where(data_flat > 0)\n            channel = np.ravel_multi_index(nonzero[:-1], shape[:-2])\n            value = data_flat[nonzero].astype(float)\n            cols = [\n                fits.Column(\"PIX\", \"J\", array=nonzero[-1]),\n                fits.Column(\"CHANNEL\", \"I\", array=channel),\n                fits.Column(\"VALUE\", \"E\", array=value),\n            ]\n        else:\n            data_flat = []\n            channel = []\n            pix = []\n            for i, _ in np.ndenumerate(npix[0]):\n                data_i = np.ravel(data[i[::-1]])\n                data_i[~np.isfinite(data_i)] = 0\n                pix_i = np.where(data_i > 0)\n                data_i = data_i[pix_i]\n                data_flat += [data_i]\n                pix += pix_i\n                channel += [\n                    np.ones(data_i.size, dtype=int)\n                    * np.ravel_multi_index(i[::-1], shape[:-2])\n                ]\n\n            pix = np.concatenate(pix)\n            channel = np.concatenate(channel)\n            value = np.concatenate(data_flat).astype(float)\n\n            cols = [\n                fits.Column(\"PIX\", \"J\", array=pix),\n                fits.Column(\"CHANNEL\", \"I\", array=channel),\n                fits.Column(\"VALUE\", \"E\", array=value),\n            ]\n\n        return fits.BinTableHDU.from_columns(cols, header=header, name=hdu)\n", "path": "gammapy/maps/wcsmap.py"}], "after_files": [{"content": "# Licensed under a 3-clause BSD style license - see LICENSE.rst\nimport json\nimport numpy as np\nfrom astropy.io import fits\nfrom .core import Map\nfrom .utils import find_bands_hdu, find_hdu, JsonQuantityEncoder\nfrom .wcs import WcsGeom\n\n\n__all__ = [\"WcsMap\"]\n\n\ndef identify_wcs_format(hdu):\n    if hdu is None:\n        return \"gadf\"\n    elif hdu.name == \"ENERGIES\":\n        return \"fgst-template\"\n    elif hdu.name == \"EBOUNDS\":\n        return \"fgst-ccube\"\n    else:\n        return \"gadf\"\n\n\nclass WcsMap(Map):\n    \"\"\"Base class for WCS map classes.\n\n    Parameters\n    ----------\n    geom : `~gammapy.maps.WcsGeom`\n        A WCS geometry object.\n    data : `~numpy.ndarray`\n        Data array.\n    \"\"\"\n\n    @classmethod\n    def create(\n        cls,\n        map_type=\"wcs\",\n        npix=None,\n        binsz=0.1,\n        width=None,\n        proj=\"CAR\",\n        frame=\"icrs\",\n        refpix=None,\n        axes=None,\n        skydir=None,\n        dtype=\"float32\",\n        meta=None,\n        unit=\"\",\n    ):\n        \"\"\"Factory method to create an empty WCS map.\n\n        Parameters\n        ----------\n        map_type : {'wcs', 'wcs-sparse'}\n            Map type.  Selects the class that will be used to\n            instantiate the map.\n        npix : int or tuple or list\n            Width of the map in pixels. A tuple will be interpreted as\n            parameters for longitude and latitude axes.  For maps with\n            non-spatial dimensions, list input can be used to define a\n            different map width in each image plane.  This option\n            supersedes width.\n        width : float or tuple or list\n            Width of the map in degrees.  A tuple will be interpreted\n            as parameters for longitude and latitude axes.  For maps\n            with non-spatial dimensions, list input can be used to\n            define a different map width in each image plane.\n        binsz : float or tuple or list\n            Map pixel size in degrees.  A tuple will be interpreted\n            as parameters for longitude and latitude axes.  For maps\n            with non-spatial dimensions, list input can be used to\n            define a different bin size in each image plane.\n        skydir : tuple or `~astropy.coordinates.SkyCoord`\n            Sky position of map center.  Can be either a SkyCoord\n            object or a tuple of longitude and latitude in deg in the\n            coordinate system of the map.\n        frame : {\"icrs\", \"galactic\"}, optional\n            Coordinate system, either Galactic (\"galactic\") or Equatorial (\"icrs\").\n        axes : list\n            List of non-spatial axes.\n        proj : string, optional\n            Any valid WCS projection type. Default is 'CAR' (cartesian).\n        refpix : tuple\n            Reference pixel of the projection.  If None then this will\n            be chosen to be center of the map.\n        dtype : str, optional\n            Data type, default is float32\n        meta : `dict`\n            Dictionary to store meta data.\n        unit : str or `~astropy.units.Unit`\n            The unit of the map\n\n        Returns\n        -------\n        map : `~WcsMap`\n            A WCS map object.\n        \"\"\"\n        from .wcsnd import WcsNDMap\n\n        geom = WcsGeom.create(\n            npix=npix,\n            binsz=binsz,\n            width=width,\n            proj=proj,\n            skydir=skydir,\n            frame=frame,\n            refpix=refpix,\n            axes=axes,\n        )\n\n        if map_type == \"wcs\":\n            return WcsNDMap(geom, dtype=dtype, meta=meta, unit=unit)\n        elif map_type == \"wcs-sparse\":\n            raise NotImplementedError\n        else:\n            raise ValueError(f\"Invalid map type: {map_type!r}\")\n\n    @classmethod\n    def from_hdulist(cls, hdu_list, hdu=None, hdu_bands=None, format=\"gadf\"):\n        \"\"\"Make a WcsMap object from a FITS HDUList.\n\n        Parameters\n        ----------\n        hdu_list :  `~astropy.io.fits.HDUList`\n            HDU list containing HDUs for map data and bands.\n        hdu : str\n            Name or index of the HDU with the map data.\n        hdu_bands : str\n            Name or index of the HDU with the BANDS table.\n        format : {'gadf', 'fgst-ccube', 'fgst-template'}\n            FITS format convention.\n\n        Returns\n        -------\n        wcs_map : `WcsMap`\n            Map object\n        \"\"\"\n        if hdu is None:\n            hdu = find_hdu(hdu_list)\n        else:\n            hdu = hdu_list[hdu]\n\n        if hdu_bands is None:\n            hdu_bands = find_bands_hdu(hdu_list, hdu)\n\n        if hdu_bands is not None:\n            hdu_bands = hdu_list[hdu_bands]\n\n        format = identify_wcs_format(hdu_bands)\n\n        wcs_map = cls.from_hdu(hdu, hdu_bands, format=format)\n\n        if wcs_map.unit.is_equivalent(\"\"):\n            if format == \"fgst-template\":\n                if \"GTI\" in hdu_list:  # exposure maps have an additional GTI hdu\n                    wcs_map.unit = \"cm2 s\"\n                else:\n                    wcs_map.unit = \"cm-2 s-1 MeV-1 sr-1\"\n\n        return wcs_map\n\n    def to_hdulist(self, hdu=None, hdu_bands=None, sparse=False, format=\"gadf\"):\n        \"\"\"Convert to `~astropy.io.fits.HDUList`.\n\n        Parameters\n        ----------\n        hdu : str\n            Name or index of the HDU with the map data.\n        hdu_bands : str\n            Name or index of the HDU with the BANDS table.\n        sparse : bool\n            Sparsify the map by only writing pixels with non-zero\n            amplitude.\n        format : {'gadf', 'fgst-ccube','fgst-template'}\n            FITS format convention.\n\n        Returns\n        -------\n        hdu_list : `~astropy.io.fits.HDUList`\n\n        \"\"\"\n        if sparse:\n            hdu = \"SKYMAP\" if hdu is None else hdu.upper()\n        else:\n            hdu = \"PRIMARY\" if hdu is None else hdu.upper()\n\n        if sparse and hdu == \"PRIMARY\":\n            raise ValueError(\"Sparse maps cannot be written to the PRIMARY HDU.\")\n\n        if format in [\"fgst-ccube\", \"fgst-template\"]:\n            if self.geom.axes[0].name != \"energy\" or len(self.geom.axes) > 1:\n                raise ValueError(\n                    \"All 'fgst' formats don't support extra axes except for energy.\"\n                )\n\n        if hdu_bands is None:\n            hdu_bands = f\"{hdu.upper()}_BANDS\"\n\n        if self.geom.axes:\n            hdu_bands_out = self.geom.to_bands_hdu(\n                hdu_bands=hdu_bands, format=format\n            )\n            hdu_bands = hdu_bands_out.name\n        else:\n            hdu_bands = None\n\n        hdu_out = self.to_hdu(hdu=hdu, hdu_bands=hdu_bands, sparse=sparse)\n\n        hdu_out.header[\"META\"] = json.dumps(self.meta, cls=JsonQuantityEncoder)\n\n        hdu_out.header[\"BUNIT\"] = self.unit.to_string(\"fits\")\n\n        if hdu == \"PRIMARY\":\n            hdulist = [hdu_out]\n        else:\n            hdulist = [fits.PrimaryHDU(), hdu_out]\n\n        if self.geom.axes:\n            hdulist += [hdu_bands_out]\n\n        return fits.HDUList(hdulist)\n\n    def to_hdu(self, hdu=\"SKYMAP\", hdu_bands=None, sparse=False):\n        \"\"\"Make a FITS HDU from this map.\n\n        Parameters\n        ----------\n        hdu : str\n            The HDU extension name.\n        hdu_bands : str\n            The HDU extension name for BANDS table.\n        sparse : bool\n            Set INDXSCHM to SPARSE and sparsify the map by only\n            writing pixels with non-zero amplitude.\n\n        Returns\n        -------\n        hdu : `~astropy.io.fits.BinTableHDU` or `~astropy.io.fits.ImageHDU`\n            HDU containing the map data.\n        \"\"\"\n        header = self.geom.to_header()\n\n        if self.is_mask:\n            data = self.data.astype(int)\n        else:\n            data = self.data\n\n        if hdu_bands is not None:\n            header[\"BANDSHDU\"] = hdu_bands\n\n        if sparse:\n            hdu_out = self._make_hdu_sparse(data, self.geom.npix, hdu, header)\n        elif hdu == \"PRIMARY\":\n            hdu_out = fits.PrimaryHDU(data, header=header)\n        else:\n            hdu_out = fits.ImageHDU(data, header=header, name=hdu)\n\n        return hdu_out\n\n    @staticmethod\n    def _make_hdu_sparse(data, npix, hdu, header):\n        shape = data.shape\n\n        # We make a copy, because below we modify `data` to handle non-finite entries\n        # TODO: The code below could probably be simplified to use expressions\n        # that create new arrays instead of in-place modifications\n        # But first: do we want / need the non-finite entry handling at all and always cast to 64-bit float?\n        data = data.copy()\n\n        if len(shape) == 2:\n            data_flat = np.ravel(data)\n            non_zero = np.where(~(data_flat == 0))\n            value = data_flat[non_zero].astype(float)\n            cols = [\n                fits.Column(\"PIX\", \"J\", array=non_zero[0]),\n                fits.Column(\"VALUE\", \"E\", array=value),\n            ]\n        elif npix[0].size == 1:\n            shape_flat = shape[:-2] + (shape[-1] * shape[-2],)\n            data_flat = np.ravel(data).reshape(shape_flat)\n            nonzero = np.where(~(data_flat == 0))\n            channel = np.ravel_multi_index(nonzero[:-1], shape[:-2])\n            value = data_flat[nonzero].astype(float)\n            cols = [\n                fits.Column(\"PIX\", \"J\", array=nonzero[-1]),\n                fits.Column(\"CHANNEL\", \"I\", array=channel),\n                fits.Column(\"VALUE\", \"E\", array=value),\n            ]\n        else:\n            data_flat = []\n            channel = []\n            pix = []\n            for i, _ in np.ndenumerate(npix[0]):\n                data_i = np.ravel(data[i[::-1]])\n                pix_i = np.where(~(data_i == 0))\n                data_i = data_i[pix_i]\n                data_flat += [data_i]\n                pix += pix_i\n                channel += [\n                    np.ones(data_i.size, dtype=int)\n                    * np.ravel_multi_index(i[::-1], shape[:-2])\n                ]\n\n            pix = np.concatenate(pix)\n            channel = np.concatenate(channel)\n            value = np.concatenate(data_flat).astype(float)\n\n            cols = [\n                fits.Column(\"PIX\", \"J\", array=pix),\n                fits.Column(\"CHANNEL\", \"I\", array=channel),\n                fits.Column(\"VALUE\", \"E\", array=value),\n            ]\n\n        return fits.BinTableHDU.from_columns(cols, header=header, name=hdu)\n", "path": "gammapy/maps/wcsmap.py"}]}