diff --git a/ChangeLog b/ChangeLog index 43eed90..584ec67 100644 --- a/ChangeLog +++ b/ChangeLog @@ -3,6 +3,7 @@ * Version 1.1.0 released ========================= + Handle FITS extensions in GSkyMap load() and save() methods (#1695) Restrict Monte Carlo simulations in GModelSpatialDiffuseConst to simulation cone (#1674) Use default GTime in GModelSpectral eval() and eval_gradients() methods (#1672) Optimise GModelSpatialDiffuseMap::set_mc_cone() method (#1673) diff --git a/NEWS b/NEWS index 98eacba..b39eabc 100644 --- a/NEWS +++ b/NEWS @@ -151,7 +151,11 @@ None 12. Sky module -------------- -None +The GSkyMap load() and save() methods now support the handling of FITS +extension names. Users can now specify with which extension name a sky +map is written to a FITS file, and from which extension name it should +be loaded. This allows to have multiple sky maps in the same file, and +to mix also HEALPix and WCS maps in the same file (#1695). 13. Support module diff --git a/include/GSkyMap.hpp b/include/GSkyMap.hpp index de28ff4..f4391d0 100644 --- a/include/GSkyMap.hpp +++ b/include/GSkyMap.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GSkyMap.hpp - Sky map class * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2015 by Juergen Knoedlseder * + * copyright (C) 2010-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -33,20 +33,12 @@ #include "GSkyDir.hpp" #include "GSkyPixel.hpp" #include "GSkyProjection.hpp" -/* -#include "GFits.hpp" -#include "GFitsTable.hpp" -#include "GFitsBinTable.hpp" -#include "GFitsImage.hpp" -#include "GFitsImageDouble.hpp" -*/ -//#include "GMatrix.hpp" -//#include "GVector.hpp" #include "GBilinear.hpp" /* __ Forward declarations _______________________________________________ */ class GFilename; class GFits; +class GFitsHDU; class GFitsTable; class GFitsBinTable; class GFitsImage; @@ -186,6 +178,8 @@ private: const GSkyDir& dir3) const; double solidangle(const GSkyDir& dir1, const GSkyDir& dir2, const GSkyDir& dir3, const GSkyDir& dir4) const; + bool is_healpix(const GFitsHDU& hdu) const; + bool is_wcs(const GFitsHDU& hdu) const; // Private data area int m_num_pixels; //!< Number of pixels (used for pixel allocation) diff --git a/src/sky/GSkyMap.cpp b/src/sky/GSkyMap.cpp index bd83ef7..9375ba8 100644 --- a/src/sky/GSkyMap.cpp +++ b/src/sky/GSkyMap.cpp @@ -1746,56 +1746,85 @@ void GSkyMap::stack_maps(void) ***************************************************************************/ void GSkyMap::load(const GFilename& filename) { - // Free memory and initialise members - free_members(); - init_members(); - - // Open FITS file - GFits fits(filename); - - // Get number of HDUs - int num = fits.size(); + // Clear sky map + clear(); // Initialize load flag - bool loaded = false; + bool loaded(false); - // First search for HEALPix extension. We can skip the first extension - // since this is always an image and a HEALPix map is stored in a - // binary table - for (int extno = 1; extno < num; ++extno) { - - // Get reference to HDU - const GFitsHDU& hdu = *fits.at(extno); + // Open FITS file + GFits fits(filename); - // If PIXTYPE keyword equals "HEALPIX" then load map - if (hdu.has_card("PIXTYPE") && hdu.string("PIXTYPE") == "HEALPIX") { + // If an extension name is specified then first try loading that + // extension + if (filename.has_extname()) { + const GFitsHDU& hdu = *fits.at(filename.extname()); + if (is_healpix(hdu)) { read_healpix(static_cast(hdu)); loaded = true; - break; } + else if (is_wcs(hdu)) { + read_wcs(static_cast(hdu)); + loaded = true; + } + } - } // endfor: looped over HDUs + // ... otherwise is an extension number is specified then try loading + // the corresponding HDU + else if (filename.has_extno()) { + const GFitsHDU& hdu = *fits.at(filename.extno()); + if (is_healpix(hdu)) { + read_healpix(static_cast(hdu)); + loaded = true; + } + else if (is_wcs(hdu)) { + read_wcs(static_cast(hdu)); + loaded = true; + } + } - // If we have not found a HEALPIX map then search now for image. - // Skip empty images + // If no map has yet been loaded then scan the file for an appropriate + // HDU and read sky map from the first suitable HDU if (!loaded) { - for (int extno = 0; extno < num; ++extno) { - // Get referene to HDU + // Get number of HDUs + int num = fits.size(); + + // First search for HEALPix extension. We can skip the first + // extension since this is always an image and a HEALPix map + // is stored in a binary table + for (int extno = 1; extno < num; ++extno) { + + // Get reference to HDU const GFitsHDU& hdu = *fits.at(extno); - // Skip if extension is not an image - if (extno > 0) { - if (hdu.string("XTENSION") != "IMAGE") - continue; + // If HDU is HEALPix then read data + if (is_healpix(hdu)) { + read_healpix(static_cast(hdu)); + loaded = true; + break; } - // Load WCS map - read_wcs(static_cast(hdu)); - break; - } // endfor: looped over HDUs - } // endif: no HEALPix map found + + // If we have not found a HEALPIX map then search now for an + // image. + if (!loaded) { + for (int extno = 0; extno < num; ++extno) { + + // Get referene to HDU + const GFitsHDU& hdu = *fits.at(extno); + + // If HDU is WCS then read data + if (is_wcs(hdu)) { + read_wcs(static_cast(hdu)); + loaded = true; + break; + } + + } // endfor: looped over HDUs + } // endif: no sky map yet loaded + } // endif: no sky map yet loaded // Close FITS file fits.close(); @@ -1806,12 +1835,14 @@ void GSkyMap::load(const GFilename& filename) /***********************************************************************//** - * @brief Save skymap into FITS file. + * @brief Save sky map into FITS file * * @param[in] filename FITS file name. - * @param[in] clobber Overwrite existing file? (true=yes) + * @param[in] clobber Overwrite existing file? (default: false) * - * The method does nothing if the skymap holds no valid WCS. + * Saves the sky map into a FITS file. If the file exists already the method + * will append the sky map to the FITS file in case that the @p clobber + * parameter is set to true. ***************************************************************************/ void GSkyMap::save(const GFilename& filename, bool clobber) const { @@ -1831,15 +1862,28 @@ void GSkyMap::save(const GFilename& filename, bool clobber) const hdu = create_wcs_hdu(); } - // Create FITS file and save it to disk + // If we have a valid HDU then save it now to the FITS file if (hdu != NULL) { - GFits fits; + + // Set extension name + if (filename.has_extname()) { + hdu->extname(filename.extname()); + } + + // Open or create FITS file (without extension name since the + // requested extension may not yet exist in the file) + GFits fits(filename.url(), true); + + // Append sky map to FITS file fits.append(*hdu); - fits.saveto(filename, clobber); - } - // Delete HDU - if (hdu != NULL) delete hdu; + // Save FITS file + fits.save(clobber); + + // Delete HDU + delete hdu; + + } // endif: HDU was valid } // endif: we had data to save @@ -2674,6 +2718,51 @@ double GSkyMap::solidangle(const GSkyDir& dir1, const GSkyDir& dir2, } +/***********************************************************************//** + * @brief Check if HDU contains HEALPix data + * + * @param[in] hdu FITS Header Data Unit. + * @return True is HDU contains HEALPix data. + ***************************************************************************/ +bool GSkyMap::is_healpix(const GFitsHDU& hdu) const +{ + // Initialise flag + bool flag(false); + + // If PIXTYPE keyword equals "HEALPIX" then signal that we have + // HEALPix data + if ((hdu.exttype() != GFitsHDU::HT_IMAGE) && + (hdu.has_card("PIXTYPE")) && + (hdu.string("PIXTYPE") == "HEALPIX")) { + flag = true; + } + + // Return flag + return (flag); +} + + +/***********************************************************************//** + * @brief Check if HDU contains WCS data + * + * @param[in] hdu FITS Header Data Unit. + * @return True is HDU contains WCS data. + ***************************************************************************/ +bool GSkyMap::is_wcs(const GFitsHDU& hdu) const +{ + // Initialise flag + bool flag(false); + + // If extension is an image thn signal that we have WCS data + if (hdu.exttype() == GFitsHDU::HT_IMAGE) { + flag = true; + } + + // Return flag + return (flag); +} + + /*========================================================================== = = = Friends = diff --git a/test/test_GSky.cpp b/test/test_GSky.cpp index ca8a5ff..82dd17c 100644 --- a/test/test_GSky.cpp +++ b/test/test_GSky.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * test_GSky.cpp - Test sky module * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2015 by Juergen Knoedlseder * + * copyright (C) 2010-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -60,6 +60,7 @@ void TestGSky::set(void){ append(static_cast(&TestGSky::test_GSkyMap_wcs_construct),"Test WCS GSkyMap constructors"); append(static_cast(&TestGSky::test_GSkyMap_wcs_io),"Test WCS GSkyMap I/O"); append(static_cast(&TestGSky::test_GSkyMap),"Test GSkyMap"); + append(static_cast(&TestGSky::test_GSkyMap_io),"Test GSkyMap I/O"); append(static_cast(&TestGSky::test_GSkyRegions_io),"Test GSkyRegions"); append(static_cast(&TestGSky::test_GSkyRegionCircle_construct),"Test GSkyRegionCircle constructors"); append(static_cast(&TestGSky::test_GSkyRegionCircle_logic),"Test GSkyRegionCircle logic"); @@ -1032,6 +1033,68 @@ void TestGSky::test_GSkyMap(void) /*************************************************************************** + * @brief GSkyMap + * + * This method tests the GSkyMap class input and output. Specifically, the + * test verifies that saving of a HEALPix and WCS map in the same FITS file + * is possible, and that upon opening the correct extension is selected. + ***************************************************************************/ +void TestGSky::test_GSkyMap_io(void) +{ + // Set filename + GFilename filename("test_skymap.fits"); + filename.remove(); + + // Create HEALPix and WCS maps and save them to FITS file + GSkyMap healpix1("GAL", 4, "RING", 1); + GSkyMap healpix2("GAL", 2, "RING", 1); + GSkyMap wcs1("CAR", "GAL", 0.0, 0.0, -1.0, 1.0, 10, 10); + GSkyMap wcs2("CAR", "GAL", 0.0, 0.0, -1.0, 1.0, 5, 5); + for (int i = 0; i < healpix1.npix(); ++i) { + healpix1(i) = double(i+1); + } + for (int i = 0; i < healpix2.npix(); ++i) { + healpix2(i) = 3.0*double(i+1); + } + for (int i = 0; i < wcs1.npix(); ++i) { + wcs1(i) = double(i+1); + } + for (int i = 0; i < wcs2.npix(); ++i) { + wcs2(i) = 3.0*double(i+1); + } + healpix1.save(filename+"[HEALPIX1]", true); + wcs1.save(filename+"[WCS1]", true); + healpix2.save(filename+"[HEALPIX2]", true); + wcs2.save(filename+"[WCS2]", true); + + // Load HEALPix map and check content + GSkyMap map1(filename+"[HEALPIX1]"); + test_value(map1.npix(), healpix1.npix(), "Check number of HEALPix pixels"); + int failures1(0); + for (int i = 0; i < map1.npix(); ++i) { + if (std::abs(map1(i)-double(i+1)) > 1.0e-6) { + failures1++; + } + } + test_value(failures1, 0, "Check number of invalid HEALPix pixel values"); + + // Load WCS map and check content + GSkyMap map2(filename+"[WCS1]"); + test_value(map2.npix(), wcs1.npix(), "Check number of WCS pixels"); + int failures2(0); + for (int i = 0; i < map2.npix(); ++i) { + if (std::abs(map2(i)-double(i+1)) > 1.0e-6) { + failures2++; + } + } + test_value(failures2, 0, "Check number of invalid WCS pixel values"); + + // Exit test + return; +} + + +/*************************************************************************** * @brief GSkyRegionCircle_construct ***************************************************************************/ void TestGSky::test_GSkyRegionCircle_construct(void) diff --git a/test/test_GSky.hpp b/test/test_GSky.hpp index bcf76d2..2a8628c 100644 --- a/test/test_GSky.hpp +++ b/test/test_GSky.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * test_GSky.hpp - Test sky module * * ----------------------------------------------------------------------- * - * copyright (C) 2012-2014 by Jean-Baptiste Cayrou * + * copyright (C) 2012-2016 by Jean-Baptiste Cayrou * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -48,6 +48,7 @@ public: void test_GSkyMap_wcs_construct(void); void test_GSkyMap_wcs_io(void); void test_GSkyMap(void); + void test_GSkyMap_io(void); void test_GSkyRegions_io(void); void test_GSkyRegionCircle_construct(void); void test_GSkyRegionCircle_logic(void);