Implement FITS extension handling in GSkyMap load() and save() methods.

Jürgen Knödlseder
1 parent c4ca5ceb
Showing 6 changed files with 208 additions and 56 deletions
ChangeLog
NEWS
include/GSkyMap.hpp
src/sky/GSkyMap.cpp
test/test_GSky.cpp
test/test_GSky.hpp
@@ -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)
@@ -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
 /***************************************************************************
  *                       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)
@@ -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<const GFitsTable&>(hdu));
             loaded = true;
-            break;
         }
+        else if (is_wcs(hdu)) {
+            read_wcs(static_cast<const GFitsImage&>(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<const GFitsTable&>(hdu));
+            loaded = true;
+        }
+        else if (is_wcs(hdu)) {
+            read_wcs(static_cast<const GFitsImage&>(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<const GFitsTable&>(hdu));
+                loaded = true;
+                break;
             }
  
-            // Load WCS map
-            read_wcs(static_cast<const GFitsImage&>(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<const GFitsImage&>(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&amp; 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&amp; 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&amp; dir1, const GSkyDir&amp; 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                                 =
 /***************************************************************************
  *                  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<pfunction>(&TestGSky::test_GSkyMap_wcs_construct),"Test WCS GSkyMap constructors");
     append(static_cast<pfunction>(&TestGSky::test_GSkyMap_wcs_io),"Test WCS GSkyMap I/O");
     append(static_cast<pfunction>(&TestGSky::test_GSkyMap),"Test GSkyMap");
+    append(static_cast<pfunction>(&TestGSky::test_GSkyMap_io),"Test GSkyMap I/O");
     append(static_cast<pfunction>(&TestGSky::test_GSkyRegions_io),"Test GSkyRegions");
     append(static_cast<pfunction>(&TestGSky::test_GSkyRegionCircle_construct),"Test GSkyRegionCircle constructors");
     append(static_cast<pfunction>(&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)
 /***************************************************************************
  *                  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);