implemented WCS loading/saving

knodlseder
1 parent 54522345
Showing 9 changed files with 675 additions and 125 deletions
include/GSkymap.hpp
include/GWcs.hpp
include/GWcsCAR.hpp
include/GWcsHPX.hpp
src/sky/GSkymap.cpp
src/sky/GWcs.cpp
src/sky/GWcsCAR.cpp
src/sky/GWcsHPX.cpp
test/test_GSky.cpp
@@ -96,7 +96,10 @@ private:
     int       xy2pix(const GSkyPixel& pix) const;
     GSkyPixel pix2xy(const int& pix) const;
     void      read_healpix(const GFitsHDU* hdu);
+    void      read_wcs(const GFitsHDU* hdu);
+    void      alloc_wcs(const GFitsHDU* hdu);
     GFitsHDU* create_healpix_hdu(void);
+    GFitsHDU* create_wcs_hdu(void);
     // Private data area
     int     m_num_pixels;   //!< Number of pixels (used for pixel allocation)
@@ -54,11 +54,11 @@ public:
     virtual GWcs& operator= (const GWcs& wcs);
     // Pure virtual methods (not implemented)
-    virtual std::string type(void) const = 0;
+    virtual void read(const GFitsHDU* hdu) = 0;
+    virtual void write(GFitsHDU* hdu) const = 0;
     // Virtual methods
-    virtual void        read(const GFitsHDU* hdu);
-    virtual void        write(GFitsHDU* hdu);
+    virtual std::string type(void) const;
     virtual std::string coordsys(void) const;
     virtual void        coordsys(const std::string& coordsys);
     virtual double      omega(const int& pix) const;
@@ -79,23 +79,32 @@ protected:
     // Typdefs
     typedef void (GWcs::*_wcspf)(GVector* coord) const;
-    // Protected methods
-    GVector wcs_dir2native(GSkyDir dir) const;
-    GSkyDir wcs_native2dir(GVector native) const;
-    void    wcs_init(const double& theta0);
-    GVector wcs_getpole(const double& theta0);
-    GMatrix wcs_get_rot(void);
-    void    dump_wcs(std::ostream& os) const;
+    // Protected methods (providing services to derived classes)
+    GVector     wcs_dir2native(GSkyDir dir) const;
+    GSkyDir     wcs_native2dir(GVector native) const;
+    void        wcs_init(const double& theta0);
+    GVector     wcs_getpole(const double& theta0);
+    GMatrix     wcs_get_rot(void);
+    void        wcs_set(const std::string& coords,
+                        const double& crval1, const double& crval2,
+                        const double& crpix1, const double& crpix2,
+                        const double& cdelt1, const double& cdelt2);
+    void        wcs_read(const GFitsHDU* hdu);
+    void        wcs_write(GFitsHDU* hdu) const;
+    std::string wcs_crval1(void) const;
+    std::string wcs_crval2(void) const;
+    void        wcs_dump(std::ostream& os) const;
     // Protected data area (astr structure)
-    int     m_coordsys;     //!< 0=celestial, 1=galactic
-    int     m_reverse;      //!< Reverse axes (1=true)
-    GVector m_crval;        //!< Value of reference point
-    GVector m_crpix;        //!< Pixel of reference point (starting from 1)
-    GVector m_cdelt;        //!< Increment at reference point (deg/pixel)
-    GVector m_npole;        //!< Native coordinate of North Pole
-    GMatrix m_cd;           //!< Astrometry parameters (2x2 matrix, deg/pixel)
-    GVector m_pv2;          //!< Projection parameters (up to 21)
+    std::string m_type;     //!< WCS type
+    int         m_coordsys; //!< 0=celestial, 1=galactic
+    int         m_reverse;  //!< Reverse axes (1=true)
+    GVector     m_crval;    //!< Value of reference point
+    GVector     m_crpix;    //!< Pixel of reference point (starting from 1)
+    GVector     m_cdelt;    //!< Increment at reference point (deg/pixel)
+    GVector     m_npole;    //!< Native coordinate of North Pole
+    GMatrix     m_cd;       //!< Astrometry parameters (2x2 matrix, deg/pixel)
+    GVector     m_pv2;      //!< Projection parameters (up to 21)
     // Derived parameters
     double  m_theta0;       //!< Native latitude of the fiducial point
@@ -54,7 +54,8 @@ public:
     GWcsCAR& operator= (const GWcsCAR& wcs);
     // Implemented pure virtual methods
-    std::string type(void) const;
+    void read(const GFitsHDU* hdu);
+    void write(GFitsHDU* hdu) const;
     // Overloaded base class methods
     double omega(const GSkyPixel& pix) const;
@@ -51,9 +51,8 @@ public:
     GWcsHPX& operator= (const GWcsHPX& wcs);
     // Implemented pure virtual methods
-    std::string type(void) const;
-    void        read(const GFitsHDU* hdu);
-    void        write(GFitsHDU* hdu) const;
+    void read(const GFitsHDU* hdu);
+    void write(GFitsHDU* hdu) const;
     // Overloaded base class methods
     double      omega(const int& pix) const;
@@ -9,7 +9,6 @@
  *   the Free Software Foundation; either version 2 of the License, or     *
  *   (at your option) any later version.                                   *
  *                                                                         *
- * ----------------------------------------------------------------------- *
  ***************************************************************************/
 /**
  * @file GSkymap.cpp
@@ -29,6 +28,7 @@
 #include "GWcsHPX.hpp"
 #include "GFits.hpp"
 #include "GFitsTableDblCol.hpp"
+#include "GFitsImageDbl.hpp"
 /* __ Method name definitions ____________________________________________ */
 #define G_CONSTRUCT_HPX     "GSkymap::GSkymap(std::string,std::string,int," \
@@ -46,6 +46,8 @@
 #define G_SET_WCS "GSkymap::set_wcs(std::string,std::string,double,double," \
                                "double,double,double,double,GMatrix,GVector)"
 #define G_READ_HEALPIX                     "GSkymap::read_healpix(GFitsHDU*)"
+#define G_READ_WCS                             "GSkymap::read_wcs(GFitsHDU*)"
+#define G_ALLOC_WCS                           "GSkymap::alloc_wcs(GFitsHDU*)"
 /* __ Macros _____________________________________________________________ */
@@ -53,6 +55,7 @@
 /* __ Debug definitions __________________________________________________ */
 //#define G_READ_HEALPIX_DEBUG                          // Debug read_healpix
+//#define G_READ_WCS_DEBUG                                  // Debug read_wcs
 /* __ Prototype __________________________________________________________ */
@@ -141,7 +144,7 @@ GSkymap::GSkymap(const std::string&amp; wcs, const std::string&amp; coords,
 /***********************************************************************//**
- * @brief 2D sky map constructor
+ * @brief WCS sky map constructor
  *
  * @param[in] wcs World Coordinate System.
  * @param[in] coords Coordinate System (CEL or GAL).
@@ -405,7 +408,7 @@ void GSkymap::load(const std::string&amp; filename)
         // Get pointer to HDU
         hdu = fits.hdu(extno);
-        // Check if PIXTYPE keyword equals "HEALPIX"
+        // If PIXTYPE keyword equals "HEALPIX" then load map
         try {
             if (hdu->card("PIXTYPE")->string() == "HEALPIX") {
                 read_healpix(hdu);
@@ -422,9 +425,23 @@ void GSkymap::load(const std::string&amp; filename)
     // Skip empty images
     if (loaded == 0) {
         for (int extno = 0; extno < num; ++extno) {
-            //TODO: Implement non HEALPix loading
-        }
-    }
+
+            // Get pointer to HDU
+            hdu = fits.hdu(extno);
+
+            // Skip if extension is not an image
+            if (extno > 0) {
+                if (hdu->card("XTENSION")->string() != "IMAGE")
+                    continue;
+            }
+
+            // Load WCS map
+            read_wcs(hdu);
+            loaded = 1;
+            break;
+
+        } // endfor: looped over HDUs
+    } // endif: no HEALPix map found
     // Close FITS file
     fits.close();
@@ -453,9 +470,8 @@ void GSkymap::save(const std::string&amp; filename, int clobber)
             hdu = create_healpix_hdu();
         // Case B: Skymap is not Healpix
-        else {
-            //TODO: Implement non HEALPix saving
-        }
+        else
+            hdu = create_wcs_hdu();
         // Create FITS file and save it to disk
         if (hdu != NULL) {
@@ -500,8 +516,14 @@ void GSkymap::read(const GFitsHDU* hdu)
         // ... otherwise try loading as non HEALPix map
         if (loaded == 0) {
-            //TODO: Implement non HEALPix loading
-        }
+
+            // Load only if HDU contains an image
+            if (hdu->card("XTENSION")->string() == "IMAGE") {
+                read_wcs(hdu);
+                loaded = 1;
+            }
+
+        } // endif
     } // endif: HDU pointer was valid
@@ -1068,6 +1090,127 @@ void GSkymap::read_healpix(const GFitsHDU* hdu)
 /***********************************************************************//**
+ * @brief Read WCS image from FITS HDU
+ *
+ * @param[in] hdu FITS HDU containing the WCS image.
+ *
+ * @exception GException::skymap_bad_image_dim
+ *            WCS image has invalid dimension (naxis=2 or 3).
+ ***************************************************************************/
+void GSkymap::read_wcs(const GFitsHDU* hdu)
+{
+    // Continue only if HDU is valid
+    if (hdu != NULL) {
+
+        // Allocate WCS
+        alloc_wcs(hdu);
+
+        // Read WCS information from FITS header
+        m_wcs->read(hdu);
+
+        // Load FITS image
+        GFitsImageDbl* image = (GFitsImageDbl*)hdu->data();
+
+        // Extract map dimension and number of maps from image
+        if (image->naxis() == 2) {
+            m_num_x    = image->naxes(0);
+            m_num_y    = image->naxes(1);
+            m_num_maps = 1;
+        }
+        else if (image->naxis() >= 3) {
+            m_num_x    = image->naxes(0);
+            m_num_y    = image->naxes(1);
+            m_num_maps = image->naxes(2);
+        }
+        else
+            throw GException::skymap_bad_image_dim(G_READ_WCS, image->naxis());
+        #if defined(G_READ_WCS_DEBUG)
+        std::cout << "m_num_x=" << m_num_x << std::endl;
+        std::cout << "m_num_y=" << m_num_y << std::endl;
+        std::cout << "m_num_maps=" << m_num_maps << std::endl;
+        #endif
+
+        // Compute number of pixels
+        m_num_pixels = m_num_x * m_num_y;
+        #if defined(G_READ_WCS_DEBUG)
+        std::cout << "m_num_pixels=" << m_num_pixels << std::endl;
+        #endif
+
+        // Allocate pixels to hold the map
+        alloc_pixels();
+
+        // Read image
+        if (image->naxis() == 2) {
+            double* ptr = m_pixels;
+            for (int iy = 0; iy < m_num_y; ++iy) {
+                for (int ix = 0; ix < m_num_x; ++ix, ++ptr)
+                    *ptr = (*image)(ix,iy);
+            }
+        }
+        else {
+            double* ptr = m_pixels;
+            for (int iy = 0; iy < m_num_y; ++iy) {
+                for (int ix = 0; ix < m_num_x; ++ix) {
+                    for (int imap = 0; imap < m_num_maps; ++imap, ++ptr)
+                         *ptr = (*image)(ix,iy,imap);
+                }
+            }
+        }
+
+    } // endif: HDU was valid
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Allocate WCS class
+ *
+ * @param[in] hdu FITS HDU containing the WCS image.
+ *
+ * @exception GException::fits_key_not_found
+ *            Unable to find required FITS header keyword.
+ * @exception GException::skymap_bad_ctype
+ *            CTYPE1 and CTYPE2 keywords are incompatible.
+ * @exception GException::wcs_invalid
+ *            WCS projection of FITS file not supported by GammaLib.
+ ***************************************************************************/
+void GSkymap::alloc_wcs(const GFitsHDU* hdu)
+{
+    // Continue only if HDU is valid
+    if (hdu != NULL) {
+
+        // Get standard keywords
+        std::string ctype1 = hdu->card("CTYPE1")->string();
+        std::string ctype2 = hdu->card("CTYPE2")->string();
+
+        // Extract projection type
+        std::string xproj = ctype1.substr(5,3);
+        std::string yproj = ctype2.substr(5,3);
+
+        // Check that projection type is identical on both axes
+        if (xproj != yproj)
+            throw GException::skymap_bad_ctype(G_ALLOC_WCS,
+                                               ctype1, ctype2);
+
+        // Allocate WCS class
+        if (xproj == "CAR") {
+            m_wcs = new GWcsCAR;
+        }
+        else
+            throw GException::wcs_invalid(G_ALLOC_WCS, xproj,
+                                          "WCS projection found in FITS file"
+                                          " is currently not supported.");
+
+    } // endif: HDU was valid
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Create FITS HDU containing Healpix data
  *
  * Returns pointer to HDU that contains the Healpix data. Deallocation of the
@@ -1123,6 +1266,71 @@ GFitsHDU* GSkymap::create_healpix_hdu(void)
 }
+/***********************************************************************//**
+ * @brief Create FITS HDU containing WCS image
+ *
+ * Returns pointer to HDU that contains the WCS image. Deallocation of the
+ * GFitsHDU object has to be done by the client.
+ ***************************************************************************/
+GFitsHDU* GSkymap::create_wcs_hdu(void)
+{
+    // Initialise result to NULL pointer
+    GFitsHDU* hdu = NULL;
+
+    // Compute size of Healpix data
+    int size = m_num_pixels * m_num_maps;
+
+    // Continue only if we have pixels
+    if (size > 0) {
+
+        // Set axis parameters for image construction
+        int naxis   = (m_num_maps == 1) ? 2 : 3;
+        int naxes[] = {m_num_x, m_num_y, m_num_maps};
+
+        // Allocate image
+        GFitsImageDbl image(naxis, naxes);
+
+        // Store data in image
+        if (naxis == 2) {
+            double* ptr = m_pixels;
+            for (int iy = 0; iy < m_num_y; ++iy) {
+                for (int ix = 0; ix < m_num_x; ++ix, ++ptr)
+                    image(ix,iy) = *ptr;
+            }
+        }
+        else {
+            double* ptr = m_pixels;
+            for (int iy = 0; iy < m_num_y; ++iy) {
+                for (int ix = 0; ix < m_num_x; ++ix) {
+                    for (int imap = 0; imap < m_num_maps; ++imap, ++ptr)
+                        image(ix,iy,imap) = *ptr;
+                }
+            }
+        }
+
+        // Create HDU that contains WCS map
+        hdu = new GFitsHDU(image);
+
+    } // endif: there were pixels
+
+    // ... otherwise create an empty header
+    else
+        hdu = new GFitsHDU;
+
+    // Set extension name
+    hdu->extname("IMAGE");
+
+    // If we have WCS information then write into FITS header
+    if (m_wcs != NULL) m_wcs->write(hdu);
+
+    // Set additional keywords
+    //TODO
+
+    // Return HDU
+    return hdu;
+}
+
+
 /*==========================================================================
  =                                                                         =
  =                              GSkymap friends                            =
@@ -25,13 +25,17 @@
 #include "GWcs.hpp"
 /* __ Method name definitions ____________________________________________ */
-#define G_WCS   "GWcs::GWcs(std::string,double,double,double,double,double,"\
-                                                                    "double)"
 #define G_COORDSYS_SET                          "GWcs::coordsys(std::string)"
 #define G_OMEGA1                                           "GWcs::omega(int)"
 #define G_OMEGA2                                     "GWcs::omega(GSkyPixel)"
 #define G_XY2DIR                                    "GWcs::xy2dir(GSkyPixel)"
 #define G_DIR2XY                                      "GWcs::dir2xy(GSkyDir)"
+#define G_WCS_SET      "GWcs::GWcs(std::string,double,double,double,double,"\
+                                                             "double,double)"
+#define G_WCS_READ                                "GWcs::wcs_read(GFitsHDU*)"
+#define G_WCS_WRITE                              "GWcs::wcs_write(GFitsHDU*)"
+#define G_WCS_CRVAL1                                     "GWcs::wcs_crval1()"
+#define G_WCS_CRVAL2                                     "GWcs::wcs_crval1()"
 /* __ Macros _____________________________________________________________ */
@@ -67,7 +71,7 @@ GWcs::GWcs(void)
 /***********************************************************************//**
- * @brief Standard 2D sky map constructor
+ * @brief Standard WCS sky map constructor
  *
  * @param[in] coords Coordinate system.
  * @param[in] crval1 X value of reference pixel.
@@ -77,7 +81,7 @@ GWcs::GWcs(void)
  * @param[in] cdelt1 Increment in x direction at reference pixel (deg).
  * @param[in] cdelt2 Increment in y direction at reference pixel (deg).
  *
- * Construct standard 2D sky map from standard definition parameters. This
+ * Construct standard WCS sky map from standard definition parameters. This
  * method
  ***************************************************************************/
 GWcs::GWcs(const std::string& coords,
@@ -89,35 +93,11 @@ GWcs::GWcs(const std::string&amp; coords,
     // Initialise class members
     init_members();
-    //TODO: Check parameters
-
-    // Set coordinate system
-    coordsys(coords);
-
-    // Set parameters
+    // Initialise reverse flag
     m_reverse = 0;
-    m_crval   = GVector(crval1, crval2);
-    m_crpix   = GVector(crpix1, crpix2);
-    m_cdelt   = GVector(cdelt1, cdelt2);
-    m_refval  = m_crval;
-    m_refpix  = m_crpix - GVector(1.0,1.0);
-    // Set CD matrix corresponding to no rotation
-    m_cd(0,0) = m_cdelt(0);
-    m_cd(1,0) = 0.0;
-    m_cd(0,1) = 0.0;
-    m_cd(1,1) = m_cdelt(1);
-
-    // Compute inverse CD matrix
-    double delta = m_cd(0,0)*m_cd(1,1) - m_cd(0,1)*m_cd(1,0);
-    if (delta != 0.0) {
-        m_invcd(0,0) =  m_cd(1,1) / delta;
-        m_invcd(0,1) = -m_cd(0,1) / delta;
-        m_invcd(1,0) = -m_cd(1,0) / delta;
-        m_invcd(1,1) =  m_cd(0,0) / delta;
-    }
-    else
-        throw GException::wcs::wcs(G_WCS, "Unable to invert CD matrix.");
+    // Set standard parameters
+    wcs_set(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
     // Return
     return;
@@ -194,37 +174,12 @@ GWcs&amp; GWcs::operator= (const GWcs&amp; wcs)
  ==========================================================================*/
 /***********************************************************************//**
- * @brief Read WCS definiton from FITS header.
- *
- * @param[in] hdu FITS HDU containing the WCS definition.
+ * @brief Return WCS type
  ***************************************************************************/
-void GWcs::read(const GFitsHDU* hdu)
+std::string GWcs::type(void) const
 {
-    // Free memory and initialise members
-    free_members();
-    init_members();
-
-    //TODO: Implement WCS definition reading
-
-    // Return
-    return;
-}
-
-
-/***********************************************************************//**
- * @brief Write WCS definiton into FITS HDU
- *
- * @param[in] hdu FITS HDU into which WCS definition will be written.
- ***************************************************************************/
-void GWcs::write(GFitsHDU* hdu)
-{
-    // Continue only if HDU is valid
-    if (hdu != NULL) {
-
-    } // endif: HDU was valid
-
-    // Return
-    return;
+    // Return Healix type
+    return m_type;
 }
@@ -753,11 +708,263 @@ GMatrix GWcs::wcs_get_rot(void)
 /***********************************************************************//**
+ * @brief Set standard WCS parameters
+ *
+ * @param[in] coords Coordinate system.
+ * @param[in] crval1 X value of reference pixel.
+ * @param[in] crval1 Y value of reference pixel.
+ * @param[in] crpix1 X index of reference pixel (first pixel is 1).
+ * @param[in] crpix2 Y index of reference pixel (first pixel is 1).
+ * @param[in] cdelt1 Increment in x direction at reference pixel (deg).
+ * @param[in] cdelt2 Increment in y direction at reference pixel (deg).
+ *
+ * @exception GException::wcs
+ *            Unable to invert CD matrix.
+ *
+ * Note that this method does not initialize the class members.
+ ***************************************************************************/
+void GWcs::wcs_set(const std::string& coords,
+                   const double& crval1, const double& crval2,
+                   const double& crpix1, const double& crpix2,
+                   const double& cdelt1, const double& cdelt2)
+
+{
+    //TODO: Check parameters
+
+    // Set coordinate system
+    coordsys(coords);
+
+    // Set parameters
+    m_crval   = GVector(crval1, crval2);
+    m_crpix   = GVector(crpix1, crpix2);
+    m_cdelt   = GVector(cdelt1, cdelt2);
+    m_refval  = m_crval;
+    m_refpix  = m_crpix - GVector(1.0,1.0);
+
+    // Set CD matrix corresponding to no rotation
+    m_cd(0,0) = m_cdelt(0);
+    m_cd(1,0) = 0.0;
+    m_cd(0,1) = 0.0;
+    m_cd(1,1) = m_cdelt(1);
+
+    // Compute inverse CD matrix
+    double delta = m_cd(0,0)*m_cd(1,1) - m_cd(0,1)*m_cd(1,0);
+    if (delta != 0.0) {
+        m_invcd(0,0) =  m_cd(1,1) / delta;
+        m_invcd(0,1) = -m_cd(0,1) / delta;
+        m_invcd(1,0) = -m_cd(1,0) / delta;
+        m_invcd(1,1) =  m_cd(0,0) / delta;
+    }
+    else
+        throw GException::wcs::wcs(G_WCS_SET, "Unable to invert CD matrix.");
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Read WCS definiton from FITS header.
+ *
+ * @param[in] hdu FITS HDU containing the WCS definition.
+ *
+ * @exception GException::fits_key_not_found
+ *            Unable to find required FITS header keyword.
+ * @exception GException::wcs_bad_coords
+ *            Coordinate system is not of valid type.
+ *
+ * Note that this method does not initialize the class members.
+ ***************************************************************************/
+void GWcs::wcs_read(const GFitsHDU* hdu)
+{
+    // Continue only if HDU is valid
+    if (hdu != NULL) {
+
+        // Get standard keywords
+        std::string ctype1 = hdu->card("CTYPE1")->string();
+        std::string ctype2 = hdu->card("CTYPE2")->string();
+        double      crval1 = hdu->card("CRVAL1")->real();
+        double      crval2 = hdu->card("CRVAL2")->real();
+        double      crpix1 = hdu->card("CRPIX1")->real();
+        double      crpix2 = hdu->card("CRPIX2")->real();
+        double      cdelt1 = hdu->card("CDELT1")->real();
+        double      cdelt2 = hdu->card("CDELT2")->real();
+
+        // Determine coordinate system
+        std::string coords;
+        std::string xcoord = ctype1.substr(0,4);
+        std::string ycoord = ctype2.substr(0,4);
+        if (xcoord == "RA--" && ycoord == "DEC-") {
+            m_reverse = 0;
+            coords    = "EQU";
+        }
+        else if (xcoord == "DEC-" && ycoord == "RA--") {
+            m_reverse = 1;
+            coords    = "EQU";
+        }
+        else if (xcoord == "GLON" && ycoord == "GLAT") {
+            m_reverse = 0;
+            coords    = "GAL";
+        }
+        else if (xcoord == "GLAT" && ycoord == "GLON") {
+            m_reverse = 1;
+            coords    = "GAL";
+        }
+        else
+            throw GException::wcs_bad_coords(G_WCS_READ, coordsys());
+
+        // Set standard parameters
+        wcs_set(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
+
+        // Get projection keywords
+        //TODO
+
+    } // endif: HDU was valid
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Write WCS definiton into FITS HDU
+ *
+ * @param[in] hdu FITS HDU into which WCS definition will be written.
+ *
+ * @exception GException::wcs_bad_coords
+ *            Coordinate system is not of valid type.
+ ***************************************************************************/
+void GWcs::wcs_write(GFitsHDU* hdu) const
+{
+    // Continue only if HDU is valid
+    if (hdu != NULL) {
+
+        // Set coordinate system dependent strings
+        std::string c_equinox;
+        std::string c_crval1;
+        std::string c_crval2;
+        std::string c_cdelt1;
+        std::string c_cdelt2;
+        std::string c_lonpole;
+        std::string c_latpole;
+        if (m_coordsys == 0) {
+            c_equinox = "Equinox for celestial coordinate system";
+            c_crval1  = "[deg] Right Ascension at reference pixel";
+            c_crval2  = "[deg] Declination at reference pixel";
+            c_cdelt1  = "[deg/pixel] Right Ascension increment at"
+                        " reference pixel";
+            c_cdelt2  = "[deg/pixel] Declination increment at"
+                        " reference pixel";
+            c_lonpole = "[deg] Native longitude of Celestial North Pole";
+            c_latpole = "[deg] Native latitude of Celestial North Pole";
+        }
+        else if (m_coordsys == 1) {
+            c_equinox = "Equinox for celestial coordinate system";
+            c_crval1 = "[deg] Galactic Longitude at reference pixel";
+            c_crval2 = "[deg] Galactic Latitude at reference pixel";
+            c_cdelt1  = "[deg/pixel] Galactic Longitude increment at"
+                        " reference pixel";
+            c_cdelt2  = "[deg/pixel] Galactic Latitude increment at"
+                        " reference pixel";
+            c_lonpole = "[deg] Native longitude of Galactic North Pole";
+            c_latpole = "[deg] Native latitude of Galactic North Pole";
+        }
+        else
+            throw GException::wcs_bad_coords(G_WCS_WRITE, coordsys());
+
+        // Set keywords
+        hdu->card("EQUINOX", 2000.0,       c_equinox);
+        hdu->card("CTYPE1",  wcs_crval1(), "Projection Type");
+        hdu->card("CTYPE2",  wcs_crval2(), "Projection Type");
+        hdu->card("CDELT1",  m_cdelt(0),   c_cdelt1);
+        hdu->card("CDELT2",  m_cdelt(1),   c_cdelt2);
+        hdu->card("CRPIX1",  m_crpix(0),
+                  "X index of reference pixel (starting from 1)");
+        hdu->card("CRPIX2",  m_crpix(1),
+                  "Y index of reference pixel (starting from 1)");
+        hdu->card("CRVAL1",  m_crval(0),   c_crval1);
+        hdu->card("CRVAL2",  m_crval(1),   c_crval2);
+        hdu->card("CROTA2",  0.0,          "[deg] Rotation Angle");
+        hdu->card("LONPOLE", m_npole(0),   c_lonpole);
+        hdu->card("LATPOLE", m_npole(1),   c_lonpole);
+        hdu->card("PV2_1",   0.0,          "Projection parameter 1");
+        hdu->card("PV2_2",   0.0,          "Projection parameter 2");
+
+    } // endif: HDU was valid
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Return crval1 string
+ *
+ * @exception GException::wcs_invalid
+ *            WCS projection not valid.
+ * @exception GException::wcs_bad_coords
+ *            Coordinate system is not of valid type.
+ ***************************************************************************/
+std::string GWcs::wcs_crval1(void) const
+{
+    // Initialise result
+    std::string crval;
+
+    // Check on correct type length
+    if (type().length() != 3)
+        throw GException::wcs_invalid(G_WCS_CRVAL1, type(),
+              "3-character type required.");
+
+    // Set coordinate system dependent value
+    if (m_coordsys == 0)
+        crval = "RA---" + type();
+    else if (m_coordsys == 1)
+        crval = "GLON-" + type();
+    else
+        throw GException::wcs_bad_coords(G_WCS_CRVAL1, coordsys());
+
+    // Return
+    return crval;
+}
+
+
+/***********************************************************************//**
+ * @brief Return crval2 string
+ *
+ * @exception GException::wcs_invalid
+ *            WCS projection not valid.
+ * @exception GException::wcs_bad_coords
+ *            Coordinate system is not of valid type.
+ ***************************************************************************/
+std::string GWcs::wcs_crval2(void) const
+{
+    // Initialise result
+    std::string crval;
+
+    // Check on correct type length
+    if (type().length() != 3)
+        throw GException::wcs_invalid(G_WCS_CRVAL2, type(),
+              "3-character type required.");
+
+    // Set coordinate system dependent value
+    if (m_coordsys == 0)
+        crval = "DEC--" + type();
+    else if (m_coordsys == 1)
+        crval = "GLAT-" + type();
+    else
+        throw GException::wcs_bad_coords(G_WCS_CRVAL1, coordsys());
+
+    // Return
+    return crval;
+}
+
+
+/***********************************************************************//**
  * @brief Dump WCS information into output stream
  *
  * @param[in] os Output stream into which WCS is dumped.
  ***************************************************************************/
-void GWcs::dump_wcs(std::ostream& os) const
+void GWcs::wcs_dump(std::ostream& os) const
 {
     // Put information into stream
     os << " Coordinate system .........: " << coordsys() << std::endl;
@@ -793,6 +1000,7 @@ void GWcs::dump_wcs(std::ostream&amp; os) const
 void GWcs::init_members(void)
 {
     // Initialise members
+    m_type.clear();
     m_coordsys = 0;
     m_reverse  = 0;
     m_crval    = GVector(2);
@@ -828,6 +1036,7 @@ void GWcs::init_members(void)
 void GWcs::copy_members(const GWcs& wcs)
 {
     // Copy attributes
+    m_type     = wcs.m_type;
     m_coordsys = wcs.m_coordsys;
     m_reverse  = wcs.m_reverse;
     m_crval    = wcs.m_crval;
@@ -188,12 +188,46 @@ GWcsCAR&amp; GWcsCAR::operator= (const GWcsCAR&amp; wcs)
  ==========================================================================*/
 /***********************************************************************//**
- * @brief Returns WCS type
+ * @brief Read WCS definiton from FITS header
+ *
+ * @param[in] hdu FITS HDU containing the Healpix definition.
  ***************************************************************************/
-std::string GWcsCAR::type(void) const
+void GWcsCAR::read(const GFitsHDU* hdu)
 {
-    // Return Healix type
-    return "CAR";
+    // Free memory and initialise members
+    free_members();
+    init_members();
+
+    // Continue only if HDU is valid
+    if (hdu != NULL) {
+
+        // Read standard WCS definition
+        wcs_read(hdu);
+
+    } // endif: HDU was valid
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Write Healpix definiton into FITS HDU
+ *
+ * @param[in] hdu FITS HDU to which the Healpix definition will be written.
+ ***************************************************************************/
+void GWcsCAR::write(GFitsHDU* hdu) const
+{
+    // Continue only if HDU is valid
+    if (hdu != NULL) {
+
+        // Write standard WCS definition
+        wcs_write(hdu);
+
+    } // endif: HDU was valid
+
+    // Return
+    return;
 }
@@ -223,6 +257,7 @@ double GWcsCAR::omega(const GSkyPixel&amp; pix) const
 void GWcsCAR::init_members(void)
 {
     // Initialise members
+    m_type = "CAR";
     // Return
     return;
@@ -315,7 +350,7 @@ std::ostream&amp; operator&lt;&lt; (std::ostream&amp; os, const GWcsCAR&amp; wcs)
     os << "=== GWcsCAR ===" << std::endl;
     // Add WCS information
-    wcs.dump_wcs(os);
+    wcs.wcs_dump(os);
     // Return output stream
     return os;
@@ -211,16 +211,6 @@ GWcsHPX&amp; GWcsHPX::operator= (const GWcsHPX&amp; wcs)
  ==========================================================================*/
 /***********************************************************************//**
- * @brief Returns WCS type
- ***************************************************************************/
-std::string GWcsHPX::type(void) const
-{
-    // Return Healix type
-    return "HPX";
-}
-
-
-/***********************************************************************//**
  * @brief Read Healpix definiton from FITS header
  *
  * @param[in] hdu FITS HDU containing the Healpix definition.
@@ -498,6 +488,7 @@ void GWcsHPX::ordering(const std::string&amp; ordering)
 void GWcsHPX::init_members(void)
 {
     // Initialise members
+    m_type        = "HPX";
     m_nside       = 0;
     m_npface      = 0;
     m_ncap        = 0;
@@ -267,8 +267,8 @@ void test_GSkymap_healpix_io(void)
     }
     catch (std::exception &e) {
         std::cout << std::endl 
-                  << "TEST ERROR: Unable to create Healpix instance from FITS file."
-                  << std::endl;
+                  << "TEST ERROR: Unable to create Healpix instance"
+                     " from FITS file." << std::endl;
         std::cout << e.what() << std::endl;
         throw;
     }
@@ -284,15 +284,15 @@ void test_GSkymap_healpix_io(void)
 /***************************************************************************
- *  Test: GSkymap_map_construct                                            *
+ *  Test: GSkymap_wcs_construct                                            *
  ***************************************************************************/
-void test_GSkymap_map_construct(void)
+void test_GSkymap_wcs_construct(void)
 {
     // Set precision
-    double eps = 1.0e-10;
-    
+    double eps = 1.0e-6;
+
     // Dump header
-    std::cout << "Test non Healpix GSkymap constructors: ";
+    std::cout << "Test WCS GSkymap constructors: ";
     // Test void constructor
     try {
@@ -323,19 +323,18 @@ void test_GSkymap_map_construct(void)
     // Test CAR projection
     try {
         GSkymap map1("CAR", "GAL", 138.817, 37.293, 0.521, 0.931, 100, 100);
-std::cout << map1 << std::endl;
         GSkyDir dir;
-        for (int l = -179; l < 179; ++l) {
-            for (int b = -89; b < 89; ++b) {
-                dir.lb_deg(double(l),double(b));        
+        for (int l = -180; l < 180; ++l) {
+            for (int b = -90; b < 90; ++b) {
+                dir.lb_deg(double(l),double(b));
                 GSkyPixel pixel    = map1.dir2xy(dir);
                 GSkyDir   dir_back = map1.xy2dir(pixel);
-                double    dl = dir.l_deg() - dir_back.l_deg();
-                double    db = dir.b_deg() - dir_back.b_deg();
-                if (fabs(dl) > eps || fabs(db) > eps) {
+                double    dist     = dir.dist_deg(dir_back);
+                if (dist > eps) {
                     std::cout << std::endl
                       << "TEST ERROR: Sky direction differs: dir="
-                      << dir << " dir_back=" << dir_back << std::endl;
+                      << dir << " dir_back=" << dir_back << " dist="
+                      << dist << " deg" << std::endl;
                     throw;
                 }
             }
@@ -361,6 +360,101 @@ std::cout &lt;&lt; map1 &lt;&lt; std::endl;
 /***************************************************************************
+ *  Test: GSkymap_wcs_io                                                   *
+ ***************************************************************************/
+void test_GSkymap_wcs_io(void)
+{
+    // Set filenames
+    const std::string file1 = "test_skymap_wcs_1.fits";
+    const std::string file2 = "test_skymap_wcs_2.fits";
+
+    // Dump header
+    std::cout << "Test WCS GSkymap I/O: ";
+
+    // Define WCS map for comparison
+    GSkymap refmap1("CAR", "GAL", 0.0, 0.0, -1.0, 1.0, 100, 100);
+    GSkymap refmap2("CAR", "CEL", 0.0, 0.0, -0.3, 0.3, 200, 200);
+
+    // Test WCS map saving
+    try {
+        for (int pix = 0; pix < refmap1.npix(); ++pix)
+            refmap1(pix) = pix+1;
+        refmap1.save(file1, 1);
+        for (int pix = 0; pix < refmap2.npix(); ++pix)
+            refmap2(pix) = pix+1;
+        refmap2.save(file2, 1);
+    }
+    catch (std::exception &e) {
+        std::cout << std::endl 
+                  << "TEST ERROR: Unable to save WCS map."
+                  << std::endl;
+        std::cout << e.what() << std::endl;
+        throw;
+    }
+    std::cout << ".";
+
+    // Test WCS map loading (1)
+    try {
+        GSkymap map;
+        map.load(file1);
+        int diff = 0;
+        for (int pix = 0; pix < refmap1.npix(); ++pix) {
+            if (map(pix) != refmap1(pix))
+                diff++;
+        }
+        if (diff > 0) {
+            std::cout << std::endl
+                      << "TEST ERROR: Loaded file "+file1+
+                         " differs from saved file ."
+                      << std::endl;
+            throw;
+        }
+    }
+    catch (std::exception &e) {
+        std::cout << std::endl 
+                  << "TEST ERROR: Unable to load WCS map."
+                  << std::endl;
+        std::cout << e.what() << std::endl;
+        throw;
+    }
+    std::cout << ".";
+
+    // Test WCS map loading (2)
+    try {
+        GSkymap map;
+        map.load(file2);
+        int diff = 0;
+        for (int pix = 0; pix < refmap2.npix(); ++pix) {
+            if (map(pix) != refmap2(pix))
+                diff++;
+        }
+        if (diff > 0) {
+            std::cout << std::endl
+                      << "TEST ERROR: Loaded file "+file2+
+                         " differs from saved file ."
+                      << std::endl;
+            throw;
+        }
+    }
+    catch (std::exception &e) {
+        std::cout << std::endl 
+                  << "TEST ERROR: Unable to load WCS map."
+                  << std::endl;
+        std::cout << e.what() << std::endl;
+        throw;
+    }
+    std::cout << ".";
+
+    // Signal final test success
+    std::cout << " ok." << std::endl;
+
+    // Exit test
+    return;
+
+}
+
+
+/***************************************************************************
  *                            Main test function                           *
  ***************************************************************************/
 int main(void)
@@ -374,7 +468,8 @@ int main(void)
     // Execute Healpix tests
     test_GSkymap_healpix_construct();
     test_GSkymap_healpix_io();
-    test_GSkymap_map_construct();
+    test_GSkymap_wcs_construct();
+    test_GSkymap_wcs_io();
     // Return
     return 0;