Fix TS map header writing and temp file handling

CHANGES.rst

@@ -100,6 +100,8 @@ Pull requests
 - Using consistent plotting style in docs [#317] (Axel Donath) 
 - Add an "About Gammapy" page to the docs [#312] (Christoph Deil)
 - Fix Debian install instructions [#326] (Victor Zabalza)
+- Fixed writing TS map headers [#273] (Axel Donath)
+- Changed temporary file handling in tests [#273] (Axel Donath)
 
 .. _gammapy_0p2_release:
 

gammapy/background/kernel.py

@@ -1,8 +1,9 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 from __future__ import print_function, division
+import os
+
 import numpy as np
 from astropy.io import fits
-from astropy.wcs import WCS
 import gc
 from ..stats import significance
 from ..image import binary_dilation_circle
@@ -196,16 +197,16 @@ def save_files(self, filebase, index):
 
         header = self.header
 
-        filename = filebase + '{0:02d}_mask'.format(index) + '.fits'
+        filename = os.path.join(filebase, '{0:02d}_mask.fits'.format(index))
         hdu = fits.ImageHDU(data=self._data[-1].mask.astype(np.uint8),
                             header=header)
         hdu.writeto(filename, clobber=True)
 
-        filename = filebase + '{0:02d}_background'.format(index) + '.fits'
+        filename = os.path.join(filebase, '{0:02d}_background.fits'.format(index))
         hdu = fits.ImageHDU(data=self._data[-1].background, header=header)
         hdu.writeto(filename, clobber=True)
 
-        filename = filebase + '{0:02d}_significance'.format(index) + '.fits'
+        filename = os.path.join(filebase, '{0:02d}_significance.fits'.format(index))
         hdu = fits.ImageHDU(data=self._data[-1].significance, header=header)
         hdu.writeto(filename, clobber=True)
 

gammapy/background/tests/test_kernel.py

@@ -1,9 +1,8 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 from __future__ import print_function, division
 import numpy as np
-from numpy.testing import assert_allclose, assert_equal
+from numpy.testing import assert_allclose
 import os
-import tempfile
 from astropy.io import fits
 from astropy.tests.helper import pytest
 from astropy.units import Quantity
@@ -12,7 +11,6 @@
 from ...image import make_empty_image
 from ...stats import significance
 from ...datasets import FermiGalacticCenter
-from ...irf import EnergyDependentTablePSF
 
 
 try:
@@ -138,7 +136,6 @@ def test_run_iteration_blob(self):
         self.ibe_blob.run_iteration()
         # Should be run twice to update the mask
         self.ibe_blob.run_iteration()
-        mask = self.ibe_blob.mask_image_hdu.data
         background = self.ibe_blob.background_image_hdu.data
         # Check background, should be 42 uniformly within 10%
         assert_allclose(background, 42 * np.ones((10, 10)), rtol=0.15)
@@ -150,16 +147,15 @@ def test_run(self):
         assert_allclose(mask.sum(), 97)
         assert_allclose(background, 42 * np.ones((10, 10)))
 
-    def test_save_files(self):
+    def test_save_files(self, tmpdir):
         """Tests that files are saves, and checks values within them."""
         # Create temporary file to write output into
-        dir = tempfile.mkdtemp()
         self.ibe.run_iteration(1)
-        self.ibe.save_files(filebase=dir, index=0)
+        self.ibe.save_files(filebase=str(tmpdir), index=0)
 
-        mask_filename = dir + '00_mask.fits'
-        significance_filename = dir + '00_significance.fits'
-        background_filename = dir + '00_background.fits'
+        mask_filename = str(tmpdir.join('00_mask.fits'))
+        significance_filename = str(tmpdir.join('00_significance.fits'))
+        background_filename = str(tmpdir.join('00_background.fits'))
 
         mask_data = fits.open(mask_filename)[1].data
         significance_data = fits.open(significance_filename)[1].data
@@ -169,5 +165,3 @@ def test_save_files(self):
         assert_allclose(mask_data.sum(), 97)
         assert_allclose(significance_data.sum(), 157.316195729298)
         assert_allclose(background_data.sum(), 4200)
-
-        os.removedirs(dir)

gammapy/background/tests/test_models.py

@@ -1,11 +1,9 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 from __future__ import print_function, division
-from tempfile import NamedTemporaryFile
 import numpy as np
 from numpy.testing import assert_allclose
 from astropy.tests.helper import pytest, remote_data, assert_quantity_allclose
 from astropy.table import Table
-from astropy.utils.data import get_pkg_data_filename
 from astropy.units import Quantity
 from astropy.coordinates import Angle
 from astropy.modeling.models import Gaussian1D
@@ -112,13 +110,13 @@ def test_spectrum_plot(self):
         assert_allclose(plot_data[:, 1], model_data.value)
 
     @remote_data
-    def test_write_fits_table(self):
+    def test_write_fits_table(self, tmpdir):
 
         filename = datasets.get_path('../test_datasets/background/bg_cube_model_test.fits',
                                      location='remote')
         bg_model_1 = CubeBackgroundModel.read(filename, format='table')
 
-        outfile = NamedTemporaryFile(suffix='.fits').name
+        outfile = str(tmpdir.join('cubebackground_table_test.fits'))
         bg_model_1.write(outfile, format='table')
 
         # test if values are correct in the saved file: compare both files
@@ -133,13 +131,13 @@ def test_write_fits_table(self):
                                  bg_model_1.energy_bins)
 
     @remote_data
-    def test_read_write_fits_image(self):
+    def test_read_write_fits_image(self, tmpdir):
 
         filename = datasets.get_path('../test_datasets/background/bg_cube_model_test.fits',
                                      location='remote')
         bg_model_1 = CubeBackgroundModel.read(filename, format='table')
 
-        outfile = NamedTemporaryFile(suffix='.fits').name
+        outfile = str(tmpdir.join('cubebackground_image_test.fits'))
         bg_model_1.write(outfile, format='image')
 
         # test if values are correct in the saved file: compare both files

gammapy/data/tests/test_spectral_cube.py

@@ -1,7 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 from __future__ import (absolute_import, division, print_function,
                         unicode_literals)
-from tempfile import NamedTemporaryFile
 import numpy as np
 from numpy.testing import assert_allclose
 from astropy.coordinates import Angle
@@ -44,13 +43,13 @@ def test_init(self):
         spectral_cube = SpectralCube(data, wcs, energy)
         assert spectral_cube.data.shape == (30, 21, 61)
 
-    def test_read_write(self):
+    def test_read_write(self, tmpdir):
         data = self.spectral_cube.data
         wcs = self.spectral_cube.wcs
         energy = self.spectral_cube.energy
 
         spectral_cube = SpectralCube(data, wcs, energy)
-        outfile = NamedTemporaryFile(suffix='.fits').name
+        outfile = str(tmpdir.join('spectral_cube_test.fits'))
         spectral_cube.writeto(outfile)
 
         spectral_cube.read(outfile)
@@ -225,8 +224,6 @@ def make_test_cubes(energies, nxpix, nypix, binsz):
     spectral_cube : `~gammapy.spectral_cube.SpectralCube`
         Cube of differential fluxes in units of cm^-2 s^-1 GeV^-1 sr^-1
     """
-    hdu = make_empty_image(nxpix, nypix, binsz)
-    solid_angle_array = solid_angle(hdu)
     header = make_header(nxpix, nypix, binsz)
     header['NAXIS'] = 3
     header['NAXIS3'] = len(energies)

gammapy/datasets/data/poisson_stats_image/simulate_data.py

@@ -2,33 +2,51 @@
 """Simulate test data (see README.md) with Sherpa"""
 from __future__ import print_function, division
 import numpy as np
-import sherpa.astro.ui as sau
+
+from astropy.modeling.models import Gaussian2D, Const2D
+from astropy.io import fits
+from astropy.wcs import WCS
+
+from gammapy.utils.random import get_random_state
 
 # Define width of the source and the PSF
 sigma_psf, sigma_source = 3, 4
-# for relation of sigma and fwhm see
-# http://cxc.harvard.edu/sherpa/ahelp/gauss2d.html
-sigma_to_fwhm = np.sqrt(8 * np.log(2))  # ~ 2.35
 sigma = np.sqrt(sigma_psf ** 2 + sigma_source ** 2)
-fwhm = sigma_to_fwhm * sigma
+amplitude = 1E3 / (2 * np.pi * sigma ** 2)
+
+source = Gaussian2D(amplitude, 99, 99, sigma, sigma)
+background = Const2D(1)
+model = source + background
+
+# Define data shape
+shape = (200, 200)
+y, x = np.indices(shape)
+
+# Create a new WCS object
+w = WCS(naxis=2)
+
+# Set up an Galactic projection
+w.wcs.crpix = [99, 99]
+w.wcs.cdelt = np.array([0.02, 0.02])
+w.wcs.crval = [0, 0]
+w.wcs.ctype = ['GLON-CAR', 'GLAT-CAR']
+
+# Fake data
+random_state = get_random_state(0)
+data = random_state.poisson(model(x, y))
+
+# Save data
+header = w.to_header()
 
-# Seed the random number generator to make the output reproducible
-np.random.seed(0)
+hdu = fits.PrimaryHDU(data=data.astype('int16'), header=header)
+hdu.writeto('counts.fits.gz', clobber=True)
 
-sau.dataspace2d((200, 200))
-sau.set_source('normgauss2d.source + const2d.background')
-sau.set_par('source.xpos', 100)
-sau.set_par('source.ypos', 100)
-sau.set_par('source.ampl', 1e3)
-sau.set_par('source.fwhm', fwhm)
-sau.set_par('background.c0', 1)
+hdu = fits.PrimaryHDU(data=model(x, y).astype('float32'), header=header)
+hdu.writeto('model.fits.gz', clobber=True)
 
-sau.fake()
-sau.save_model('model.fits.gz', clobber=True)
-sau.save_data('counts.fits.gz', clobber=True)
+hdu = fits.PrimaryHDU(data=background(x, y).astype('int16'), header=header)
+hdu.writeto('background.fits.gz', clobber=True)
 
-sau.set_source('source')
-sau.save_model('source.fits.gz', clobber=True)
+hdu = fits.PrimaryHDU(data=source(x, y).astype('float32'), header=header)
+hdu.writeto('source.fits.gz', clobber=True)
 
-sau.set_source('background')
-sau.save_model('background.fits.gz', clobber=True)

gammapy/datasets/load.py

@@ -162,14 +162,15 @@ def load_poisson_stats_image(extra_info=False, return_filenames=False):
                 out[name] = filename
             else:
                 data = fits.getdata(filename)
-                out[name] = data
+                out[name] = data.astype('float64')
     else:
         filename = get_path('poisson_stats_image/counts.fits.gz')
         if return_filenames:
             out = filename
         else:
-            out = fits.getdata(filename)
-
+            out = fits.getdata(filename).astype('float64')
+    if extra_info and not return_filenames:
+        out['header'] = fits.getheader(get_path('poisson_stats_image/counts.fits.gz'))
     return out
 
 

gammapy/datasets/tests/test_catalogs.py

@@ -1,44 +1,43 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 from __future__ import (absolute_import, division, print_function,
                         unicode_literals)
-from tempfile import NamedTemporaryFile
 from astropy.tests.helper import remote_data
 from ... import datasets
 
 
 @remote_data
-def test_load_catalog_atnf():
+def test_load_catalog_atnf(tmpdir):
     catalog = datasets.load_catalog_atnf()
     assert len(catalog) == 2399
 
     # Check if catalog can be serialised to FITS
-    filename = NamedTemporaryFile(suffix='.fits').name
+    filename = str(tmpdir.join('atnf_test.fits'))
     catalog.write(filename)
 
 
 # TODO: activate test when available
 @remote_data
-def _test_load_catalog_hess_galactic():
+def _test_load_catalog_hess_galactic(tmpdir):
     catalog = datasets.load_catalog_hess_galactic()
     assert len(catalog) == 42
 
     # Check if catalog can be serialised to FITS
-    filename = NamedTemporaryFile(suffix='.fits').name
+    filename = str(tmpdir.join('hess_test.fits'))
     catalog.write(filename)
 
 
 @remote_data
-def test_load_catalog_green():
+def test_load_catalog_green(tmpdir):
     catalog = datasets.load_catalog_green()
     assert len(catalog) == 294
 
     # Check if catalog can be serialised to FITS
-    filename = NamedTemporaryFile(suffix='.fits').name
+    filename = str(tmpdir.join('green_test.fits'))
     catalog.write(filename)
 
 
 @remote_data
-def test_load_catalog_snrcat():
+def test_load_catalog_snrcat(tmpdir):
     snrcat = datasets.fetch_catalog_snrcat()
 
     # Check SNR table
@@ -47,7 +46,7 @@ def test_load_catalog_snrcat():
     expected_colnames = ['Source_Name', 'RAJ2000']
     assert set(expected_colnames).issubset(table.colnames)
     # Check if catalog can be serialised to FITS
-    filename = NamedTemporaryFile(suffix='.fits').name
+    filename = str(tmpdir.join('snrcat_test.fits'))
     table.write(filename)
 
     # Check OBS table
@@ -57,15 +56,15 @@ def test_load_catalog_snrcat():
     expected_colnames = ['SNR_id', 'source_id']
     assert set(expected_colnames).issubset(table.colnames)
     # Check if catalog can be serialised to FITS
-    filename = NamedTemporaryFile(suffix='.fits').name
+    filename = str(tmpdir.join('obs_test.fits'))
     table.write(filename)
 
 
 @remote_data
-def test_load_catalog_tevcat():
+def test_load_catalog_tevcat(tmpdir):
     catalog = datasets.load_catalog_tevcat()
     assert len(catalog) == 173
 
     # Check if catalog can be serialised to FITS
-    filename = NamedTemporaryFile(suffix='.fits').name
+    filename = str(tmpdir.join('tevcat_test.fits'))
     catalog.write(filename)

gammapy/detect/test_statistics.py

@@ -67,33 +67,34 @@ def read(cls, filename):
         Read TS map result from file.
         """
         hdu_list = fits.open(filename)
-        ts = hdu_list['ts'].data
-        sqrt_ts = hdu_list['sqrt_ts'].data
-        amplitude = hdu_list['amplitude'].data
-        niter = hdu_list['niter'].data
+        ts = hdu_list['ts'].data.astype('float64')
+        sqrt_ts = hdu_list['sqrt_ts'].data.astype('float64')
+        amplitude = hdu_list['amplitude'].data.astype('float64')
+        niter = hdu_list['niter'].data.astype('float64')
         scale = hdu_list[0].header['SCALE']
         if scale == 'max':
             scale = hdu_list['scale'].data
-        morphology = hdu_list[0].header['MORPH']
+        morphology = hdu_list[0].header.get('MORPH')
         return cls(ts=ts, sqrt_ts=sqrt_ts, amplitude=amplitude, niter=niter,
                    scale=scale, morphology=morphology)
 
     def write(self, filename, header, overwrite=False):
         """Write TS map results to file"""
+        header = header.copy()
         hdu_list = fits.HDUList()
         if 'MORPH' not in header and hasattr(self, 'morphology'):
             header['MORPH'] = self.morphology, 'Source morphology assumption.'
         if not np.isscalar(self.scale):
             header['EXTNAME'] = 'scale'
             header['HDUNAME'] = 'scale'
             header['SCALE'] = 'max',  'Source morphology scale parameter.'
-            hdu_list.append(fits.ImageHDU(self.scale.astype('float32'), header))
+            hdu_list.append(fits.ImageHDU(self.scale.astype('float64'), header))
         else:
             header['SCALE'] = self.scale,  'Source morphology scale parameter.'
         for key in ['ts', 'sqrt_ts', 'amplitude', 'niter']:
             header['EXTNAME'] = key
             header['HDUNAME'] = key
-            hdu_list.append(fits.ImageHDU(self[key].astype('float32'), header))
+            hdu_list.append(fits.ImageHDU(self[key].astype('float64'), header))
 
         hdu_list.writeto(filename, clobber=overwrite)