Implement generic GSkyPixel class and rename pixel/index/direction conversion methods consistently.

The GSkyPixel class now handles 1D and 2D sky indices. This made the distinctive 1D and 2D methods in GSkyProjection and the derived classes obsolete. Now, all projection classes take GSkyPixel as argument. The GSkymap class now uses the universal GSkyPixel class. Methods have been renamed consistently. An integer index is now consistently called "sky map index" and the GSkyPixel object is named a "sky map pixel". Methods operating on indices contain now "inx", method operating on sky map pixels "pix".

Implement generic GSkyPixel class and rename pixel/index/direction conversion methods consistently.
The GSkyPixel class now handles 1D and 2D sky indices. This made the distinctive 1D and 2D methods in GSkyProjection and the derived classes obsolete. Now, all projection classes take GSkyPixel as argument. The GSkymap class now uses the universal GSkyPixel class. Methods have been renamed consistently. An integer index is now consistently called "sky map index" and the GSkyPixel object is named a "sky map pixel". Methods operating on indices contain now "inx", method operating on sky map pixels "pix".
Jürgen Knödlseder
1 parent 257bd8bf
Showing 22 changed files with 779 additions and 764 deletions
include/GHealpix.hpp
include/GSkyProjection.hpp
include/GSkymap.hpp
include/GWcs.hpp
inst/com/src/GCOMEventCube.cpp
inst/com/src/GCOMSupport.cpp
inst/cta/src/GCTAEventCube.cpp
inst/lat/src/GLATEventCube.cpp
inst/lat/src/GLATMeanPsf.cpp
inst/lat/src/GLATResponse.cpp
inst/lat/test/test_LAT.cpp
pyext/GHealpix.i
pyext/GSkyProjection.i
pyext/GSkymap.i
pyext/GWcs.i
src/model/GModelSpatialDiffuseCube.cpp
src/model/GModelSpatialDiffuseMap.cpp
src/sky/GHealpix.cpp
src/sky/GSkymap.cpp
src/sky/GWcs.cpp
@@ -61,16 +61,14 @@ public:
     // Implemented pure virtual base class methods
     virtual void        clear(void);
     virtual GHealpix*   clone(void) const;
+    virtual int         size(void) const;
     virtual std::string code(void) const;
     virtual std::string name(void) const;
     virtual void        read(const GFitsHDU* hdu);
     virtual void        write(GFitsHDU* hdu) const;
-    virtual double      omega(const int& pix) const;
-    virtual double      omega(const GSkyPixel& pix) const;
-    virtual GSkyDir     pix2dir(const int& pix) const;
-    virtual int         dir2pix(const GSkyDir& dir) const;
-    virtual GSkyDir     xy2dir(const GSkyPixel& pix) const;
-    virtual GSkyPixel   dir2xy(const GSkyDir& dir) const;
+    virtual double      omega(const GSkyPixel& pixel) const;
+    virtual GSkyDir     pix2dir(const GSkyPixel& pixel) const;
+    virtual GSkyPixel   dir2pix(const GSkyDir& dir) const;
     virtual std::string print(const GChatter& chatter = NORMAL) const;
     // Other methods
@@ -108,6 +106,20 @@ private:
 /***********************************************************************//**
+ * @brief Return dimension of projection
+ *
+ * @return Dimension of projection.
+ *
+ * Returns the dimension of the projection.
+ ***************************************************************************/
+inline
+int GHealpix::size(void) const
+{
+    return 1;
+}
+
+
+/***********************************************************************//**
  * @brief Return projection code
  *
  * @return Projection code.
@@ -62,16 +62,14 @@ public:
     // Pure virtual methods
     virtual void            clear(void) = 0;
     virtual GSkyProjection* clone(void) const = 0;
+    virtual int             size(void) const = 0;
     virtual std::string     code(void) const = 0;
     virtual std::string     name(void) const = 0;
     virtual void            read(const GFitsHDU* hdu) = 0;
     virtual void            write(GFitsHDU* hdu) const = 0;
-    virtual double          omega(const int& pix) const = 0;
-    virtual double          omega(const GSkyPixel& pix) const = 0;
-    virtual GSkyDir         pix2dir(const int& pix) const = 0;
-    virtual int             dir2pix(const GSkyDir& dir) const = 0;
-    virtual GSkyDir         xy2dir(const GSkyPixel& pix) const = 0;
-    virtual GSkyPixel       dir2xy(const GSkyDir& dir) const = 0;
+    virtual double          omega(const GSkyPixel& pixel) const = 0;
+    virtual GSkyDir         pix2dir(const GSkyPixel& pixel) const = 0;
+    virtual GSkyPixel       dir2pix(const GSkyDir& dir) const = 0;
     virtual std::string     print(const GChatter& chatter = NORMAL) const = 0;
     // Virtual methods
@@ -45,20 +45,22 @@
 /***********************************************************************//**
  * @class GSkymap
  *
- * @brief GSkymap class interface defintion
+ * @brief Sky map class
  *
- * This class implements a set of skymaps. Skymaps may either be present in
+ * This class implements a set of sky maps. Skymaps may either be present in
  * the HEALPix format or in any (supported) WCS format. Skymap pixels are
  * stored in a double precision array indexed by (x,y,map), with the x
  * axis being the most rapidely varying axis. 
  * 
- * Skymap pixels may be either accessed via 1D operators (i,map), where i 
- * runs over the (x,y) direction of the map, or via 2D operators (pixel,map),
+ * Skymap pixels may be either accessed via index operators (i,map), where i 
+ * runs over the (x,y) direction of the map, or via pixel operators (pixel,map),
  * where pixel implements a 2D (x,y) pixel direction. The first operator is
  * the preferred access method for HEALPix maps while the second operator is
  * the preferred access method for WCS maps. Conversion methods between the
  * index or sky pixel and the true physical sky direction are provided by the
  * pix2dir() and dir2pix() methods.
+ *
+ * @todo Rewrite class description.
  ***************************************************************************/
 class GSkymap : public GBase {
@@ -83,44 +85,38 @@ public:
     virtual ~GSkymap(void);
     // Operators
-    GSkymap& operator=(const GSkymap& map);
-
-    // 1D pixel methods
-    double&       operator()(const int& pixel, const int& map = 0);
-    const double& operator()(const int& pixel, const int& map = 0) const;
-    GSkyDir       pix2dir(const int& pix) const;
-    int           dir2pix(const GSkyDir& dir) const;
-    double        omega(const int& pix) const;
-
-    // 2D pixel methods
+    GSkymap&      operator=(const GSkymap& map);
+    double&       operator()(const int& index, const int& map = 0);
+    const double& operator()(const int& index, const int& map = 0) const;
     double&       operator()(const GSkyPixel& pixel, const int& map = 0);
     const double& operator()(const GSkyPixel& pixel, const int& map = 0) const;
-    GSkyDir       xy2dir(const GSkyPixel& pix) const;
-    GSkyPixel     dir2xy(const GSkyDir& dir) const;
-    double        omega(const GSkyPixel& pix) const;
-
-    // Sky direction methods
     double        operator()(const GSkyDir& dir, const int& map = 0) const;
     // Methods
-    void            clear(void);
-    GSkymap*        clone(void) const;
-    void            load(const std::string& filename);
-    void            save(const std::string& filename, bool clobber = false) const;
-    void            read(const GFitsHDU* hdu);
-    void            write(GFits* file) const;
-    int             npix(void) const;
-    int             nx(void) const;
-    int             ny(void) const;
-    int             nmaps(void) const;
-    int             xy2pix(const GSkyPixel& pix) const;
-    GSkyPixel       pix2xy(const int& pix) const;
-    GSkyProjection* projection(void) const { return m_proj; }
-    void            projection(const GSkyProjection& proj);
-    double*         pixels(void) const { return m_pixels; }
-    bool            contains(const GSkyDir& dir) const;
-    bool            contains(const GSkyPixel& pixel) const;
-    std::string     print(const GChatter& chatter = NORMAL) const;
+    void                  clear(void);
+    GSkymap*              clone(void) const;
+    const int&            npix(void) const;
+    const int&            nx(void) const;
+    const int&            ny(void) const;
+    const int&            nmaps(void) const;
+    GSkyPixel             inx2pix(const int& index) const;
+    GSkyDir               inx2dir(const int& index) const;
+    GSkyDir               pix2dir(const GSkyPixel& pixel) const;
+    int                   pix2inx(const GSkyPixel& pixel) const;
+    int                   dir2inx(const GSkyDir& dir) const;
+    GSkyPixel             dir2pix(const GSkyDir& dir) const;
+    double                omega(const int& index) const;
+    double                omega(const GSkyPixel& pixel) const;
+    bool                  contains(const GSkyDir& dir) const;
+    bool                  contains(const GSkyPixel& pixel) const;
+    const GSkyProjection* projection(void) const;
+    void                  projection(const GSkyProjection& proj);
+    const double*         pixels(void) const;
+    void                  load(const std::string& filename);
+    void                  save(const std::string& filename, bool clobber = false) const;
+    void                  read(const GFitsHDU* hdu);
+    void                  write(GFits* file) const;
+    std::string           print(const GChatter& chatter = NORMAL) const;
 private:
     // Private methods
@@ -148,4 +144,86 @@ private:
     double*         m_pixels;     //!< Pointer to skymap pixels
 };
+
+/***********************************************************************//**
+ * @brief Returns number of pixels
+ *
+ * @return Number of pixels in one sky map.
+ *
+ * Returns the number of pixels in one sky map.
+ ***************************************************************************/
+inline
+const int& GSkymap::npix(void) const
+{
+    return m_num_pixels;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns number of pixels in x coordinate
+ *
+ * @return Number of pixels in the X direction.
+ *
+ * Returns the number of pixels in the X direction. If the sky map is a
+ * one dimensional array (which is the case for the Healpix projection),
+ * the method returns 0.
+ ***************************************************************************/
+inline
+const int& GSkymap::nx(void) const
+{
+    return m_num_x;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns number of pixels in y coordinate
+ *
+ * @return Number of pixels in the Y direction.
+ *
+ * Returns the number of pixels in the Y direction. If the sky map is a
+ * one dimensional array (which is the case for the Healpix projection),
+ * the method returns 0.
+ ***************************************************************************/
+inline
+const int& GSkymap::ny(void) const
+{
+    return m_num_y;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns number of maps
+ *
+ * @return Number of maps in the sky map object.
+ ***************************************************************************/
+inline
+const int& GSkymap::nmaps(void) const
+{
+    return m_num_maps;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns pointer to sky projection
+ *
+ * @return Pointer to sky projection (NULL if no projection is defined).
+ ***************************************************************************/
+inline
+const GSkyProjection* GSkymap::projection(void) const
+{
+    return m_proj;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns pointer to pixel data
+ *
+ * @return Pointer to pixel data.
+ ***************************************************************************/
+inline
+const double* GSkymap::pixels(void) const
+{
+    return m_pixels;
+}
+
 #endif /* GSKYMAP_HPP */
@@ -67,14 +67,12 @@ public:
     virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
     // Implemented virtual base class methods
+    virtual int         size(void) const;
     virtual void        read(const GFitsHDU* hdu);
     virtual void        write(GFitsHDU* hdu) const;
-    virtual double      omega(const int& pix) const;
-    virtual double      omega(const GSkyPixel& pix) const;
-    virtual GSkyDir     pix2dir(const int& pix) const;
-    virtual int         dir2pix(const GSkyDir& dir) const;
-    virtual GSkyDir     xy2dir(const GSkyPixel& pix) const;
-    virtual GSkyPixel   dir2xy(const GSkyDir& dir) const;
+    virtual double      omega(const GSkyPixel& pixel) const;
+    virtual GSkyDir     pix2dir(const GSkyPixel& pixel) const;
+    virtual GSkyPixel   dir2pix(const GSkyDir& dir) const;
     // Other methods
     void   set(const std::string& coords,
@@ -208,4 +206,18 @@ protected:
     //struct spcprm spc
 };
+
+/***********************************************************************//**
+ * @brief Return dimension of projection
+ *
+ * @return Dimension of projection.
+ *
+ * Returns the dimension of the projection.
+ ***************************************************************************/
+inline
+int GWcs::size(void) const
+{
+    return 2;
+}
+
 #endif /* GWCS_HPP */
@@ -504,7 +504,7 @@ int GCOMEventCube::number(void) const
     double number = 0.0;
     // Get pointer on skymap pixels
-    double* pixels = m_map.pixels();
+    const double* pixels = m_map.pixels();
     // Sum event cube
     if (size() > 0 && pixels != NULL) {
@@ -708,7 +708,7 @@ void GCOMEventCube::set_scatter_directions(void)
     for (int iy = 0; iy < npsi(); ++iy) {
         for (int ix = 0; ix < nchi(); ++ix) {
             GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
-            m_dirs.push_back(m_map.xy2dir(pixel));
+            m_dirs.push_back(m_map.pix2dir(pixel));
             m_omega.push_back(m_map.omega(pixel));
         }
     }
@@ -871,7 +871,7 @@ void GCOMEventCube::set_bin(const int&amp; index)
     m_dir.phibar(m_phi[iphi]);
     // Set pointers
-    m_bin.m_counts = &(m_map.pixels()[index]);
+    m_bin.m_counts = const_cast<double*>(&(m_map.pixels()[index]));
     m_bin.m_omega  = &(m_omega[ipix]);
     // Return
@@ -54,7 +54,7 @@
 void com_wcs_mer2car(GSkymap& map)
 {
     // Get WCS poiunter
-    GWcs* wcs = dynamic_cast<GWcs*>(map.projection());
+    const GWcs* wcs = dynamic_cast<const GWcs*>(map.projection());
     // Apply kluge only if the skymap is in Mercator's projection
     if (wcs != NULL && wcs->code() == "MER") {
@@ -418,7 +418,7 @@ int GCTAEventCube::number(void) const
     double number = 0.0;
     // Get pointer on skymap pixels
-    double* pixels = m_map.pixels();
+    const double* pixels = m_map.pixels();
     // Sum event cube
     if (size() > 0 && pixels != NULL) {
@@ -670,7 +670,7 @@ void GCTAEventCube::set_directions(void)
     for (int iy = 0; iy < ny(); ++iy) {
         for (int ix = 0; ix < nx(); ++ix) {
             GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
-            m_dirs.push_back(GCTAInstDir(m_map.xy2dir(pixel)));
+            m_dirs.push_back(GCTAInstDir(m_map.pix2dir(pixel)));
             m_omega.push_back(m_map.omega(pixel));
         }
     }
@@ -770,24 +770,27 @@ void GCTAEventCube::set_bin(const int&amp; index)
 {
     // Optionally check if the index is valid
     #if defined(G_RANGE_CHECK)
-    if (index < 0 || index >= size())
+    if (index < 0 || index >= size()) {
         throw GException::out_of_range(G_SET_BIN, index, 0, size()-1);
+    }
     #endif
     // Check for the existence of energies and energy widths
-    if (m_energies.size() != ebins() || m_ewidth.size() != ebins())
+    if (m_energies.size() != ebins() || m_ewidth.size() != ebins()) {
         throw GCTAException::no_energies(G_SET_BIN);
+    }
     // Check for the existence of sky directions and solid angles
-    if (m_dirs.size() != npix() || m_omega.size() != npix())
+    if (m_dirs.size() != npix() || m_omega.size() != npix()) {
         throw GCTAException::no_dirs(G_SET_BIN);
+    }
     // Get pixel and energy bin indices.
     int ipix = index % npix();
     int ieng = index / npix();
     // Set pointers
-    m_bin.m_counts = &(m_map.pixels()[index]);
+    m_bin.m_counts = const_cast<double*>(&(m_map.pixels()[index]));
     m_bin.m_energy = &(m_energies[ieng]);
     m_bin.m_time   = &m_time;
     m_bin.m_dir    = &(m_dirs[ipix]);
@@ -388,7 +388,7 @@ int GLATEventCube::number(void) const
     double number = 0.0;
     // Get pointer on skymap pixels
-    double* pixels = m_map.pixels();
+    const double* pixels = m_map.pixels();
     // Sum event cube
     if (size() > 0 && pixels != NULL) {
@@ -561,7 +561,7 @@ double GLATEventCube::maxrad(const GSkyDir&amp; srcDir) const
     int iy = 0;
     for (int ix = 0; ix < nx(); ++ix) {
         GSkyPixel pixel    = GSkyPixel(double(ix), double(iy));
-        double    distance = m_map.xy2dir(pixel).dist_deg(srcDir);
+        double    distance = m_map.pix2dir(pixel).dist_deg(srcDir);
         if (distance < radius) {
             radius = distance;
         }
@@ -571,7 +571,7 @@ double GLATEventCube::maxrad(const GSkyDir&amp; srcDir) const
     iy = ny()-1;
     for (int ix = 0; ix < nx(); ++ix) {
         GSkyPixel pixel    = GSkyPixel(double(ix), double(iy));
-        double    distance = m_map.xy2dir(pixel).dist_deg(srcDir);
+        double    distance = m_map.pix2dir(pixel).dist_deg(srcDir);
         if (distance < radius) {
             radius = distance;
         }
@@ -581,7 +581,7 @@ double GLATEventCube::maxrad(const GSkyDir&amp; srcDir) const
     int ix = 0;
     for (int iy = 0; iy < ny(); ++iy) {
         GSkyPixel pixel    = GSkyPixel(double(ix), double(iy));
-        double    distance = m_map.xy2dir(pixel).dist_deg(srcDir);
+        double    distance = m_map.pix2dir(pixel).dist_deg(srcDir);
         if (distance < radius) {
             radius = distance;
         }
@@ -591,7 +591,7 @@ double GLATEventCube::maxrad(const GSkyDir&amp; srcDir) const
     ix = nx()-1;
     for (int iy = 0; iy < ny(); ++iy) {
         GSkyPixel pixel    = GSkyPixel(double(ix), double(iy));
-        double    distance = m_map.xy2dir(pixel).dist_deg(srcDir);
+        double    distance = m_map.pix2dir(pixel).dist_deg(srcDir);
         if (distance < radius) {
             radius = distance;
         }
@@ -841,7 +841,7 @@ void GLATEventCube::set_directions(void)
     for (int iy = 0; iy < ny(); ++iy) {
         for (int ix = 0; ix < nx(); ++ix) {
             GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
-            m_dirs.push_back(GLATInstDir(m_map.xy2dir(pixel)));
+            m_dirs.push_back(GLATInstDir(m_map.pix2dir(pixel)));
             m_omega.push_back(m_map.omega(pixel));
         }
     }
@@ -967,7 +967,7 @@ void GLATEventCube::set_bin(const int&amp; index)
     // Set pointers
     m_bin.m_cube   = this;
-    m_bin.m_counts = &(m_map.pixels()[index]);
+    m_bin.m_counts = const_cast<double*>(&(m_map.pixels()[index]));
     m_bin.m_energy = &(m_energies[m_bin.m_ieng]);
     m_bin.m_time   = &m_time;
     m_bin.m_dir    = &(m_dirs[m_bin.m_ipix]);
@@ -738,7 +738,7 @@ void GLATMeanPsf::set_map_corrections(const GLATObservation&amp; obs)
                     // Compute offset angle in degrees
                     GSkyPixel pixel  = GSkyPixel(double(ix), double(iy));
-                    double    offset = cube->map().xy2dir(pixel).dist_deg(m_dir);
+                    double    offset = cube->map().pix2dir(pixel).dist_deg(m_dir);
                     double    omega  = cube->map().omega(pixel);
                     // Use only pixels within maximum PSF radius
@@ -400,10 +400,10 @@ double GLATResponse::irf(const GLATEventBin&amp; event,
         nodes.set_value(srcEng.log10MeV());
         // Compute diffuse response
-        GSkymap* map    = cube->diffrsp(idiff);
-        double*  pixels = map->pixels() + event.ipix();
-        rsp             = nodes.wgt_left()  * pixels[nodes.inx_left()  * map->npix()] +
-                          nodes.wgt_right() * pixels[nodes.inx_right() * map->npix()];
+        GSkymap*      map    = cube->diffrsp(idiff);
+        const double* pixels = map->pixels() + event.ipix();
+        rsp                  = nodes.wgt_left()  * pixels[nodes.inx_left()  * map->npix()] +
+                               nodes.wgt_right() * pixels[nodes.inx_right() * map->npix()];
         // Divide by solid angle and ontime since source maps are given in units of
         // counts/pixel/MeV.
@@ -287,7 +287,7 @@ void TestGLATLtCube::test_one_ltcube(const std::string&amp; datadir, const double&amp; r
         GLATLtCube ltcube(lat_ltcube);
         GEnergy energy;
         for (int i = 0; i < map.npix(); ++i) {
-            GSkyDir dir = map.pix2dir(i);
+            GSkyDir dir = map.inx2dir(i);
             map(i) = ltcube(dir, energy, test_fct1);
         }
         map.save(file1, true);
@@ -304,7 +304,7 @@ void TestGLATLtCube::test_one_ltcube(const std::string&amp; datadir, const double&amp; r
         GLATLtCube ltcube(lat_ltcube);
         GEnergy energy;
         for (int i = 0; i < map.npix(); ++i) {
-            GSkyDir dir = map.pix2dir(i);
+            GSkyDir dir = map.inx2dir(i);
             map(i) = ltcube(dir, energy, test_fct2);
         }
         map.save(file2, true);
@@ -47,16 +47,14 @@ public:
     // Implemented pure virtual methods
     virtual void        clear(void);
     virtual GHealpix*   clone(void) const;
+    virtual int         size(void) const;
     virtual std::string code(void) const;
     virtual std::string name(void) const;
     virtual void        read(const GFitsHDU* hdu);
     virtual void        write(GFitsHDU* hdu) const;
-    virtual double      omega(const int& pix) const;
-    virtual double      omega(const GSkyPixel& pix) const;
-    virtual GSkyDir     pix2dir(const int& pix) const;
-    virtual int         dir2pix(const GSkyDir& dir) const;
-    virtual GSkyDir     xy2dir(const GSkyPixel& pix) const;
-    virtual GSkyPixel   dir2xy(const GSkyDir& dir) const;
+    virtual double      omega(const GSkyPixel& pixel) const;
+    virtual GSkyDir     pix2dir(const GSkyPixel& pixel) const;
+    virtual GSkyPixel   dir2pix(const GSkyDir& dir) const;
     // Other methods
     const int&   npix(void) const;
@@ -48,16 +48,14 @@ public:
     // Pure virtual methods (not implemented)
     virtual void            clear(void) = 0;
     virtual GSkyProjection* clone(void) const = 0;
+    virtual int             size(void) const = 0;
     virtual std::string     code(void) const = 0;
     virtual std::string     name(void) const = 0;
     virtual void            read(const GFitsHDU* hdu) = 0;
     virtual void            write(GFitsHDU* hdu) const = 0;
-    virtual double          omega(const int& pix) const = 0;
-    virtual double          omega(const GSkyPixel& pix) const = 0;
-    virtual GSkyDir         pix2dir(const int& pix) const = 0;
-    virtual int             dir2pix(const GSkyDir& dir) const = 0;
-    virtual GSkyDir         xy2dir(const GSkyPixel& pix) const = 0;
-    virtual GSkyPixel       dir2xy(const GSkyDir& dir) const = 0;
+    virtual double          omega(const GSkyPixel& pixel) const = 0;
+    virtual GSkyDir         pix2dir(const GSkyPixel& pixel) const = 0;
+    virtual GSkyPixel       dir2pix(const GSkyDir& dir) const = 0;
     // Virtual methods
     virtual std::string coordsys(void) const;
@@ -127,34 +127,30 @@ public:
     GSkymap(const GSkymap& map);
     virtual ~GSkymap(void);
-    // 1D pixel methods
-    GSkyDir       pix2dir(const int& pix) const;
-    int           dir2pix(const GSkyDir& dir) const;
-    double        omega(const int& pix) const;
-
-    // 2D pixel methods
-    GSkyDir       xy2dir(const GSkyPixel& pix) const;
-    GSkyPixel     dir2xy(const GSkyDir& dir) const;
-    double        omega(const GSkyPixel& pix) const;
-
     // Methods
-    void            clear(void);
-    GSkymap*        clone(void) const;
-    void            load(const std::string& filename);
-    void            save(const std::string& filename, bool clobber = false) const;
-    void            read(const GFitsHDU* hdu);
-    void            write(GFits* file) const;
-    int             npix(void) const;
-    int             nx(void) const;
-    int             ny(void) const;
-    int             nmaps(void) const;
-    int             xy2pix(const GSkyPixel& pix) const;
-    GSkyPixel       pix2xy(const int& pix) const;
-    GSkyProjection* projection(void) const;
-    void            projection(const GSkyProjection& proj);
-    double*         pixels(void) const;
-    bool            contains(const GSkyDir& dir) const;
-    bool            contains(const GSkyPixel& pixel) const;
+    void                  clear(void);
+    GSkymap*              clone(void) const;
+    const int&            npix(void) const;
+    const int&            nx(void) const;
+    const int&            ny(void) const;
+    const int&            nmaps(void) const;
+    GSkyPixel             inx2pix(const int& index) const;
+    GSkyDir               inx2dir(const int& index) const;
+    GSkyDir               pix2dir(const GSkyPixel& pixel) const;
+    int                   pix2inx(const GSkyPixel& pixel) const;
+    int                   dir2inx(const GSkyDir& dir) const;
+    GSkyPixel             dir2pix(const GSkyDir& dir) const;
+    double                omega(const int& index) const;
+    double                omega(const GSkyPixel& pixel) const;
+    bool                  contains(const GSkyDir& dir) const;
+    bool                  contains(const GSkyPixel& pixel) const;
+    const GSkyProjection* projection(void) const;
+    void                  projection(const GSkyProjection& proj);
+    const double*         pixels(void) const;
+    void                  load(const std::string& filename);
+    void                  save(const std::string& filename, bool clobber = false) const;
+    void                  read(const GFitsHDU* hdu);
+    void                  write(GFits* file) const;
 };
@@ -53,14 +53,12 @@ public:
     virtual std::string name(void) const = 0;
     // Implemented virtual methods
+    virtual int         size(void) const;
     virtual void        read(const GFitsHDU* hdu);
     virtual void        write(GFitsHDU* hdu) const;
-    virtual double      omega(const int& pix) const;
-    virtual double      omega(const GSkyPixel& pix) const;
-    virtual GSkyDir     pix2dir(const int& pix) const;
-    virtual int         dir2pix(const GSkyDir& dir) const;
-    virtual GSkyDir     xy2dir(const GSkyPixel& pix) const;
-    virtual GSkyPixel   dir2xy(const GSkyDir& dir) const;
+    virtual double      omega(const GSkyPixel& pixel) const;
+    virtual GSkyDir     pix2dir(const GSkyPixel& pixel) const;
+    virtual GSkyPixel   dir2pix(const GSkyDir& dir) const;
     // Other methods
     void   set(const std::string& coords,
@@ -446,15 +446,15 @@ GSkyDir GModelSpatialDiffuseCube::mc(const GEnergy&amp; energy,
             }
         }
-        // Convert 1D pixel index to 2D pixel index
-        GSkyPixel pixel = m_cube.pix2xy(low-offset);
+        // Convert sky map index to sky map pixel
+        GSkyPixel pixel = m_cube.inx2pix(low-offset);
         // Randomize pixel
         pixel.x(pixel.x() + ran.uniform() - 0.5);
         pixel.y(pixel.y() + ran.uniform() - 0.5);
         // Get sky direction
-        dir = m_cube.xy2dir(pixel);
+        dir = m_cube.pix2dir(pixel);
     } // endif: there were pixels in sky map
@@ -350,15 +350,15 @@ GSkyDir GModelSpatialDiffuseMap::mc(const GEnergy&amp; energy,
             }
         }
-        // Convert 1D pixel index to 2D pixel index
-        GSkyPixel pixel = m_map.pix2xy(low);
+        // Convert sky map index to sky map pixel
+        GSkyPixel pixel = m_map.inx2pix(low);
         // Randomize pixel
         pixel.x(pixel.x() + ran.uniform() - 0.5);
         pixel.y(pixel.y() + ran.uniform() - 0.5);
         // Get sky direction
-        dir = m_map.xy2dir(pixel);
+        dir = m_map.pix2dir(pixel);
     } // endif: there were pixels in sky map
@@ -350,27 +350,14 @@ void GHealpix::write(GFitsHDU* hdu) const
 /***********************************************************************//**
  * @brief Returns solid angle of pixel
  *
- * @param[in] pix Pixel number (0,1,...,num_pixels).
+ * @param[in] pixel Sky pixel.
+ * @return Solid angle of pixel.
  *
- * HEALPix pixels have all the same solid angle, hence the pix argument is
- * not used.
+ * Returns the solid angle of the specified @p pixel. Note that HEALPix
+ * pixels have all the same solid angle, hence the @p pixel argument is in
+ * fact not used by the method.
  ***************************************************************************/
-double GHealpix::omega(const int& pix) const
-{
-    // Return solid angle
-    return m_omega;
-}
-
-
-/***********************************************************************//**
- * @brief Returns solid angle of pixel
- *
- * @param[in] pix Sky pixel.
- *
- * HEALPix pixels have all the same solid angle, hence the pix argument is
- * not used.
- ***************************************************************************/
-double GHealpix::omega(const GSkyPixel& pix) const
+double GHealpix::omega(const GSkyPixel& pixel) const
 {
     // Return solid angle
     return m_omega;
@@ -380,28 +367,35 @@ double GHealpix::omega(const GSkyPixel&amp; pix) const
 /***********************************************************************//**
  * @brief Returns sky direction of pixel
  *
- * @param[in] pix Pixel number (0,1,...,m_num_pixels).
+ * @param[in] pixel Sky map pixel.
  ***************************************************************************/
-GSkyDir GHealpix::pix2dir(const int& pix) const
+GSkyDir GHealpix::pix2dir(const GSkyPixel& pixel) const
 {
-    // Declare result
-    GSkyDir result;
-    double  theta = 0.0;
-    double  phi   = 0.0;
+    // Throw an exception if sky map pixel is not 1D
+    if (!pixel.is1D()) {
+        std::string msg = "Sky map pixel "+pixel.print()+" is not"
+                          " 1-dimensional.\n"
+                          "Only 1-dimensional pixels are supported by the"
+                          " Healpix projection.";
+        throw GException::invalid_argument(G_XY2DIR, msg);
+    }
     // Perform ordering dependent conversion
+    double theta = 0.0;
+    double phi   = 0.0;
     switch (m_ordering) {
     case 0:
-        pix2ang_ring(pix, &theta, &phi);
+        pix2ang_ring(int(pixel), &theta, &phi);
         break;
     case 1:
-        pix2ang_nest(pix, &theta, &phi);
+        pix2ang_nest(int(pixel), &theta, &phi);
         break;
     default:
         break;
     }
     // Store coordinate system dependent result
+    GSkyDir result;
     switch (m_coordsys) {
     case 0:
         result.radec(phi, gammalib::pihalf - theta);
@@ -412,25 +406,23 @@ GSkyDir GHealpix::pix2dir(const int&amp; pix) const
     default:
         break;
     }
-
-    // Return result
+    
+    // Return
     return result;
 }
 /***********************************************************************//**
- * @brief Returns pixel for a given sky direction
+ * @brief Returns sky map pixel of sky coordinate
  *
- * @param[in] dir Sky direction.
+ * @param[in] dir Sky coordinate.
+ * @return Sky map pixel.
  ***************************************************************************/
-int GHealpix::dir2pix(const GSkyDir& dir) const
+GSkyPixel GHealpix::dir2pix(const GSkyDir& dir) const
 {
-    // Declare result
-    int    pix;
+    // Compute coordinate system dependent (z,phi)
     double z;
     double phi;
-
-    // Compute coordinate system dependent (z,phi)
     switch (m_coordsys) {
     case 0:
         z   = cos(gammalib::pihalf - dir.dec());
@@ -445,75 +437,20 @@ int GHealpix::dir2pix(const GSkyDir&amp; dir) const
     }
     // Perform ordering dependent conversion
+    int index;
     switch (m_ordering) {
     case 0:
-        pix = ang2pix_z_phi_ring(z, phi);
+        index = ang2pix_z_phi_ring(z, phi);
         break;
     case 1:
-        pix = ang2pix_z_phi_nest(z, phi);
+        index = ang2pix_z_phi_nest(z, phi);
         break;
     default:
         break;
     }
-    // Return pixel index
-    return pix;
-}
-
-
-/***********************************************************************//**
- * @brief Returns sky direction of pixel
- *
- * @param[in] pix Sky pixel.
- *
- * This dummy method throws an error when called as sky pixels are not
- * implemented for a HealPix grid. Maybe we should implement them?
- *
- * @todo Think about implementation of sky pixel for HealPix data. We may
- *       use only the first argument, and even carry a usage flag in
- *       GSkyPixel
- ***************************************************************************/
-GSkyDir GHealpix::xy2dir(const GSkyPixel& pix) const
-{
-    // Set error message
-    std::string message = "Method not defined for HPX projection.";
-
-    // Throw error
-    throw GException::wcs(G_XY2DIR, message);
-
-    // Set dummy return value
-    GSkyDir result;
-    
-    // Return
-    return result;
-}
-
-
-/***********************************************************************//**
- * @brief Returns pixel of sky direction
- *
- * @param[in] dir Sky direction.
- *
- * This dummy method throws an error when called as sky pixels are not
- * implemented for a HealPix grid. Maybe we should implement them?
- *
- * @todo Think about implementation of sky pixel for HealPix data. We may
- *       use only the first argument, and even carry a usage flag in
- *       GSkyPixel
- ***************************************************************************/
-GSkyPixel GHealpix::dir2xy(const GSkyDir& dir) const
-{
-    // Set error message
-    std::string message = "Method not defined for HPX projection.";
-
-    // Throw error
-    throw GException::wcs(G_DIR2XY2, message);
-
-    // Set dummy return value
-    GSkyPixel result;
-    
-    // Return
-    return result;
+    // Return sky map pixel
+    return (GSkyPixel(index));
 }
@@ -46,13 +46,13 @@
 #define G_OP_ACCESS_1D                        "GSkymap::operator(int&, int&)"
 #define G_OP_ACCESS_2D                  "GSkymap::operator(GSkyPixel&, int&)"
 #define G_OP_VALUE                        "GSkymap::operator(GSkyDir&, int&)"
-#define G_READ                               "GSkymap::read(const GFitsHDU*)"
-#define G_PIX2DIR                                    "GSkymap::pix2dir(int&)"
+#define G_INX2DIR                                    "GSkymap::inx2dir(int&)"
+#define G_PIX2DIR                              "GSkymap::pix2dir(GSkyPixel&)"
+#define G_DIR2INX                                "GSkymap::dir2inx(GSkyDir&)"
 #define G_DIR2PIX                                "GSkymap::dir2pix(GSkyDir&)"
-#define G_XY2DIR                                "GSkymap::xy2dir(GSkyPixel&)"
-#define G_DIR2XY                                  "GSkymap::dir2xy(GSkyDir&)"
 #define G_OMEGA1                                       "GSkymap::omega(int&)"
 #define G_OMEGA2                                 "GSkymap::omega(GSkyPixel&)"
+#define G_READ                               "GSkymap::read(const GFitsHDU*)"
 #define G_SET_WCS     "GSkymap::set_wcs(std::string&, std::string&, double&,"\
                               " double&, double&, double&, double&, double&,"\
                                                        " GMatrix&, GVector&)"
@@ -260,72 +260,104 @@ GSkymap::~GSkymap(void)
  ==========================================================================*/
 /***********************************************************************//**
- * @brief 1D pixel access operator
+ * @brief Assignment operator
+ *
+ * @param[in] map Sky map.
+ * @return Sky map.
+ ***************************************************************************/
+GSkymap& GSkymap::operator=(const GSkymap& map)
+{
+    // Execute only if object is not identical
+    if (this != &map) {
+
+        // Free members
+        free_members();
+
+        // Initialise private members for clean destruction
+        init_members();
+
+        // Copy members
+        copy_members(map);
+
+    } // endif: object was not identical
+
+    // Return this object
+    return *this;
+}
+
+
+/***********************************************************************//**
+ * @brief Pixel index access operator
  *
- * @param[in] pixel Pixel index (0,1,...,m_num_pixels).
- * @param[in] map Map index (0,1,...,m_num_maps).
+ * @param[in] index Pixel index [0,...,npix()-1].
+ * @param[in] map Map index [0,...,nmaps()-1].
+ * @return Sky map pixel value.
  *
  * @exception GException::out_of_range
  *            Pixel index and/or map index are outside valid range.
  *
  * Access sky map pixel by its index, where the most quickly varying axis is
  * the x axis of the map.
- *
- * @todo Implement proper skymap exception (actual is for matrix elements)
  ***************************************************************************/
-double& GSkymap::operator() (const int& pixel, const int& map)
+double& GSkymap::operator()(const int& index, const int& map)
 {
     // Throw an error if pixel index or map index is not in valid range
     #if defined(G_RANGE_CHECK)
-    if (pixel < 0 || pixel >= m_num_pixels) {
-        throw GException::out_of_range(G_OP_ACCESS_1D, pixel, 0, m_num_pixels-1);
+    if (index < 0 || index >= m_num_pixels) {
+        throw GException::out_of_range(G_OP_ACCESS_1D,
+                                       "Sky map pixel index",
+                                       index, m_num_pixels);
     }
     if (map < 0 || map >= m_num_maps) {
-        throw GException::out_of_range(G_OP_ACCESS_1D, map, 0, m_num_maps-1);
+        throw GException::out_of_range(G_OP_ACCESS_1D,
+                                       "Sky map map index",
+                                       map, m_num_maps);
     }
     #endif
     // Return reference to pixel value
-    return m_pixels[pixel+m_num_pixels*map];
+    return m_pixels[index+m_num_pixels*map];
 }
 /***********************************************************************//**
- * @brief 1D pixel access operator (const variant)
+ * @brief Pixel index access operator (const variant)
  *
- * @param[in] pixel Pixel index (0,1,...,m_num_pixels).
- * @param[in] map Map index (0,1,...,m_num_maps).
+ * @param[in] index Pixel index [0,...,npix()-1].
+ * @param[in] map Map index [0,...,nmaps()-1].
  *
  * @exception GException::out_of_range
  *            Pixel index and/or map index are outside valid range.
  *
  * Access sky map pixel by its index, where the most quickly varying axis is
  * the x axis of the map.
- *
- * @todo Implement proper skymap exception (actual is for matrix elements)
  ***************************************************************************/
-const double& GSkymap::operator() (const int& pixel, const int& map) const
+const double& GSkymap::operator()(const int& index, const int& map) const
 {
     // Throw an error if pixel index or map index is not in valid range
     #if defined(G_RANGE_CHECK)
-    if (pixel < 0 || pixel >= m_num_pixels) {
-        throw GException::out_of_range(G_OP_ACCESS_1D, pixel, 0, m_num_pixels-1);
+    if (index < 0 || index >= m_num_pixels) {
+        throw GException::out_of_range(G_OP_ACCESS_1D,
+                                       "Sky map pixel index",
+                                       index, m_num_pixels);
     }
     if (map < 0 || map >= m_num_maps) {
-        throw GException::out_of_range(G_OP_ACCESS_1D, map, 0, m_num_maps-1);
+        throw GException::out_of_range(G_OP_ACCESS_1D,
+                                       "Sky map map index",
+                                       map, m_num_maps);
     }
     #endif
     // Return reference to pixel value
-    return m_pixels[pixel+m_num_pixels*map];
+    return m_pixels[index+m_num_pixels*map];
 }
 /***********************************************************************//**
- * @brief 2D pixel access operator
+ * @brief Sky map pixel access operator
  *
- * @param[in] pixel Sky pixel.
- * @param[in] map Map index (0,1,...,m_num_maps).
+ * @param[in] pixel Sky map pixel.
+ * @param[in] map Map index [0,...,nmaps()-1].
  *
  * @exception GException::out_of_range
  *            Sky pixel and/or map index are outside valid range.
@@ -335,7 +367,7 @@ const double&amp; GSkymap::operator() (const int&amp; pixel, const int&amp; map) const
  *
  * @todo Implement proper skymap exception (actual is for matrix elements)
  ***************************************************************************/
-double& GSkymap::operator() (const GSkyPixel& pixel, const int& map)
+double& GSkymap::operator()(const GSkyPixel& pixel, const int& map)
 {
     // Throw an error if pixel index or map index is not in valid range
     #if defined(G_RANGE_CHECK)
@@ -345,12 +377,14 @@ double&amp; GSkymap::operator() (const GSkyPixel&amp; pixel, const int&amp; map)
                                        m_num_x-1, m_num_y-1);
     }
     if (map < 0 || map >= m_num_maps) {
-        throw GException::out_of_range(G_OP_ACCESS_2D, map, 0, m_num_maps-1);
+        throw GException::out_of_range(G_OP_ACCESS_2D,
+                                       "Sky map map index",
+                                       map, m_num_maps);
     }
     #endif
     // Get pixel index
-    int index = xy2pix(pixel);
+    int index = pix2inx(pixel);
     // Return reference to pixel value
     return m_pixels[index+m_num_pixels*map];
@@ -358,10 +392,10 @@ double&amp; GSkymap::operator() (const GSkyPixel&amp; pixel, const int&amp; map)
 /***********************************************************************//**
- * @brief 2D pixel access operator
+ * @brief Sky map pixel access operator
  *
- * @param[in] pixel Sky pixel.
- * @param[in] map Map index (0,1,...,m_num_maps).
+ * @param[in] pixel Sky map pixel.
+ * @param[in] map Map index [0,...,nmaps()-1].
  *
  * @exception GException::out_of_range
  *            Sky pixel and/or map index are outside valid range.
@@ -381,12 +415,14 @@ const double&amp; GSkymap::operator()(const GSkyPixel&amp; pixel, const int&amp; map) const
                                        m_num_x-1, m_num_y-1);
     }
     if (map < 0 || map >= m_num_maps) {
-        throw GException::out_of_range(G_OP_ACCESS_2D, map, 0, m_num_maps-1);
+        throw GException::out_of_range(G_OP_ACCESS_2D,
+                                       "Sky map map index",
+                                       map, m_num_maps);
     }
     #endif
     // Get pixel index
-    int index = xy2pix(pixel);
+    int index = pix2inx(pixel);
     // Return reference to pixel value
     return m_pixels[index+m_num_pixels*map];
@@ -397,7 +433,7 @@ const double&amp; GSkymap::operator()(const GSkyPixel&amp; pixel, const int&amp; map) const
  * @brief Return interpolated skymap value for sky direction
  *
  * @param[in] dir Sky direction.
- * @param[in] map Map index (0,1,...,m_num_maps).
+ * @param[in] map Map index [0,...,nmaps()-1].
  *
  * @exception GException::out_of_range
  *            Map index lies outside valid range.
@@ -406,12 +442,14 @@ const double&amp; GSkymap::operator()(const GSkyPixel&amp; pixel, const int&amp; map) const
  * interpolation of the neighbouring pixels. If the sky direction falls
  * outside the area covered by the skymap, a value of 0 is returned.
  ***************************************************************************/
-double GSkymap::operator() (const GSkyDir& dir, const int& map) const
+double GSkymap::operator()(const GSkyDir& dir, const int& map) const
 {
     // Throw an error if the map index is not in valid range
     #if defined(G_RANGE_CHECK)
     if (map < 0 || map >= m_num_maps) {
-        throw GException::out_of_range(G_OP_VALUE, map, 0, m_num_maps-1);
+        throw GException::out_of_range(G_OP_VALUE,
+                                       "Sky map map index",
+                                       map, m_num_maps);
     }
     #endif
@@ -419,7 +457,7 @@ double GSkymap::operator() (const GSkyDir&amp; dir, const int&amp; map) const
     double intensity = 0.0;
     // Determine sky pixel
-    GSkyPixel pixel = dir2xy(dir);
+    GSkyPixel pixel = dir2pix(dir);
     // Continue only if pixel is within the map
     if (contains(pixel)) {
@@ -476,32 +514,6 @@ double GSkymap::operator() (const GSkyDir&amp; dir, const int&amp; map) const
 }
-/***********************************************************************//**
- * @brief Assignment operator
- *
- * @param[in] map Sky map.
- ***************************************************************************/
-GSkymap& GSkymap::operator= (const GSkymap& map)
-{
-    // Execute only if object is not identical
-    if (this != &map) {
-
-        // Free members
-        free_members();
-
-        // Initialise private members for clean destruction
-        init_members();
-
-        // Copy members
-        copy_members(map);
-
-    } // endif: object was not identical
-
-    // Return this object
-    return *this;
-}
-
-
 /*==========================================================================
  =                                                                         =
  =                             Public methods                              =
@@ -511,7 +523,7 @@ GSkymap&amp; GSkymap::operator= (const GSkymap&amp; map)
  /***********************************************************************//**
  * @brief Clear instance.
  *
- * This method properly resets the object to an initial state.
+ * Resets the sky map to the initial state.
  ***************************************************************************/
 void GSkymap::clear(void)
 {
@@ -527,7 +539,9 @@ void GSkymap::clear(void)
 /***********************************************************************//**
- * @brief Clone instance
+ * @brief Clone sky map
+ *
+ * @return Pointer to deep copy of sky map.
  ***************************************************************************/
 GSkymap* GSkymap::clone(void) const
 {
@@ -535,6 +549,368 @@ GSkymap* GSkymap::clone(void) const
 }
+/***********************************************************************//**
+ * @brief Converts pixel index into sky map pixel
+ *
+ * @param[in] index Pixel index [0,...,npix()-1].
+ * @return Sky map pixel.
+ *
+ * Converts the pixel @p index into a sky map pixel (GSkyPixel).The dimension
+ * of GSkyPixel will be identical to the dimension of the sky map (i.e. a 1D
+ * sky map leads to a 1D GSkyPixel object, a 2D sky map leads to a 2D
+ * GSkyPixel object.
+ ***************************************************************************/
+GSkyPixel GSkymap::inx2pix(const int& index) const
+{
+    // Initialise sky map pixel
+    GSkyPixel pixel;
+
+    // Get x and y indices
+    double x;
+    double y;
+    if (m_num_x != 0) { //!< 2D sky map
+        pixel.x(double(index % m_num_x));
+        pixel.y(double(index / m_num_x));
+    }
+    else {              //!< 1D sky map
+        pixel.index(index);
+    }
+
+    // Return pixel
+    return pixel;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns sky direction of pixel
+ *
+ * @param[in] index Pixel index [0,...,npix()-1].
+ * @return Sky direction.
+ *
+ * @exception GException::invalid_value
+ *            No valid sky projection found.
+ *
+ * Returns sky direction for a given pixel index.
+ ***************************************************************************/
+GSkyDir GSkymap::inx2dir(const int& index) const
+{
+    // Throw error if sky projection is not valid
+    if (m_proj == NULL) {
+        std::string msg = "Sky projection has not been defined.";
+        throw GException::invalid_value(G_INX2DIR, msg);
+    }
+
+    // Determine sky direction from pixel index.
+    GSkyDir dir = m_proj->pix2dir(inx2pix(index));
+
+    // Return sky direction
+    return dir;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns sky direction of pixel
+ *
+ * @param[in] pixel Sky map pixel.
+ * @return Sky direction.
+ *
+ * @exception GException::invalid_value
+ *            No valid sky projection found.
+ * @exception GException::invalid_argument
+ *            2D sky map pixel used to access 1D projection.
+ *
+ * Returns sky direction for a given sky map @p pixel.
+ ***************************************************************************/
+GSkyDir GSkymap::pix2dir(const GSkyPixel& pixel) const
+{
+    // Throw error if WCS is not valid
+    if (m_proj == NULL) {
+        std::string msg = "Sky projection has not been defined.";
+        throw GException::invalid_value(G_PIX2DIR, msg);
+    }
+
+    // Initialise sky direction
+    GSkyDir dir;
+
+    // If pixel size matches the projection size then perform a straight
+    // forward conversion
+    if (m_proj->size() == pixel.size()) {
+        dir = m_proj->pix2dir(pixel);
+    }
+
+    // ... otherwise, if we have a 2D projection but a 1D pixel then
+    // interpret the pixel as the linear index in the pixel array
+    else if (m_proj->size() == 2) {
+        dir = m_proj->pix2dir(GSkyPixel(inx2pix(int(pixel))));
+    }
+
+    // ... otherwise we have a 1D projection but a 2D pixel. There is
+    // no unambiguous way to handle this case, hence we throw an exception
+    else {
+        std::string msg = "A 2-dimensional sky map pixel "+pixel.print()+
+                          " is used to determine the sky direction for"
+                          " the 1-dimensional sky projection \""+
+                          m_proj->name()+"\"\n"
+                          "Please specify a 1-dimensional sky map pixel.";
+        throw GException::invalid_argument(G_PIX2DIR, msg);
+    }
+
+    // Return sky direction
+    return dir;
+}
+
+
+/***********************************************************************//**
+ * @brief Converts sky map pixel into pixel index
+ *
+ * @param[in] pixel Sky map pixel.
+ * @return Pixel index [0,...,npix()-1].
+ *
+ * Converts a sky map @p pixel into the pixel index.
+ ***************************************************************************/
+int GSkymap::pix2inx(const GSkyPixel& pixel) const
+{
+    // Initialise pixel index
+    int index = 0;
+
+    // Handle 1D sky map pixel
+    if (pixel.is1D()) {
+        index = int(pixel);
+    }
+
+    // Handle 2D sky map pixel
+    else if (pixel.is2D()) {
+
+        // Get x and y indices by rounding the (x,y) values
+        int ix = int(pixel.x()+0.5);
+        int iy = int(pixel.y()+0.5);
+
+        // Set index
+        index = ix + iy * m_num_x;
+
+    }
+
+    // Return index
+    return index;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns pixel index for a given sky direction
+ *
+ * @param[in] dir Sky direction.
+ * @return Pixel index [0,...,npix()-1].
+ *
+ * @exception GException::invalid_value
+ *            No valid sky projection found.
+ *
+ * Returns sky map pixel index for a given sky direction.
+ ***************************************************************************/
+int GSkymap::dir2inx(const GSkyDir& dir) const
+{
+    // Throw error if WCS is not valid
+    if (m_proj == NULL) {
+        std::string msg = "Sky projection has not been defined.";
+        throw GException::invalid_value(G_DIR2INX, msg);
+    }
+
+    // Determine pixel index for a given sky direction
+    int index = pix2inx(m_proj->dir2pix(dir));
+
+    // Return pixel index
+    return index;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns sky map pixel for a given sky direction
+ *
+ * @param[in] dir Sky direction.
+ * @return Sky map pixel.
+ *
+ * @exception GException::invalid_value
+ *            No valid sky projection found.
+ *
+ * Returns sky map pixel for a given sky direction.
+ ***************************************************************************/
+GSkyPixel GSkymap::dir2pix(const GSkyDir& dir) const
+{
+    // Throw error if WCS is not valid
+    if (m_proj == NULL) {
+        std::string msg = "Sky projection has not been defined.";
+        throw GException::invalid_value(G_DIR2PIX, msg);
+    }
+
+    // Determine pixel for a given sky direction
+    GSkyPixel pixel = m_proj->dir2pix(dir);
+
+    // Return pixel
+    return pixel;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns solid angle of pixel
+ *
+ * @param[in] index Pixel index [0,...,npix()-1].
+ * @return Solid angle (steradians)
+ *
+ * @exception GException::invalid_value
+ *            No valid sky projection found.
+ *
+ * Returns the solid angle of the pixel with the specified @p index.
+ ***************************************************************************/
+double GSkymap::omega(const int& index) const
+{
+    // Throw error if WCS is not valid
+    if (m_proj == NULL) {
+        std::string msg = "Sky projection has not been defined.";
+        throw GException::invalid_value(G_OMEGA1, msg);
+    }
+
+    // Determine solid angle from pixel index.
+    double omega = m_proj->omega(inx2pix(index));
+
+    // Return solid angle
+    return omega;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns solid angle of pixel
+ *
+ * @param[in] pixel Sky map pixel.
+ * @return Solid angle (steradians)
+ *
+ * @exception GException::invalid_value
+ *            No valid sky projection found.
+ *
+ * Returns the solid angle of the specified sky map @p pixel.
+ ***************************************************************************/
+double GSkymap::omega(const GSkyPixel& pixel) const
+{
+    // Throw error if WCS is not valid
+    if (m_proj == NULL) {
+        std::string msg = "Sky projection has not been defined.";
+        throw GException::invalid_value(G_OMEGA2, msg);
+    }
+
+    // Initialise solid angle
+    double omega = 0.0;
+
+    // If pixel size matches the projection size then perform a straight
+    // forward solid angle determination
+    if (m_proj->size() == pixel.size()) {
+        omega = m_proj->omega(pixel);
+    }
+
+    // ... otherwise, if we have a 2D projection but a 1D pixel then
+    // interpret the pixel as the linear index in the pixel array
+    else if (m_proj->size() == 2) {
+        omega = m_proj->omega(GSkyPixel(inx2pix(int(pixel))));
+    }
+
+    // ... otherwise we have a 1D projection but a 2D pixel. There is
+    // no unambiguous way to handle this case, hence we throw an exception
+    else {
+        std::string msg = "A 2-dimensional sky map pixel "+pixel.print()+
+                          " is used to determine the solid angle for"
+                          " the 1-dimensional sky projection \""+
+                          m_proj->name()+"\"\n"
+                          "Please specify a 1-dimensional sky map pixel.";
+        throw GException::invalid_argument(G_OMEGA2, msg);
+    }
+
+    // Return solid angle
+    return omega;
+}
+
+
+/***********************************************************************//**
+ * @brief Set sky projection
+ *
+ * @param[in] proj Sky projection.
+ *
+ * Sets the projection from celestial to pixel coordinates. The method
+ * performs a deep copy of @p proj, allowing to destroy the argument after
+ * using the method.
+ *
+ * Warning: this method may corrupt the GSkymap object as it allows assigning
+ * for example a 1D projection to a 2D skymap. Please use this method only
+ * when you know what you're doing.
+ *
+ * @todo We may restrict this method to not allow changing the projection
+ * dimension.
+ ***************************************************************************/
+void GSkymap::projection(const GSkyProjection& proj)
+{
+    // Free any existing WCS
+    if (m_proj != NULL) delete m_proj;
+
+    // Clone input WCS
+    m_proj = proj.clone();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Verifies if sky direction falls in map
+ *
+ * @param[in] dir Sky direction.
+ *
+ * This method checks if the specified sky direction falls within the pixels
+ * covered by the skymap. The method uses the dir2xy method to convert the
+ * sky direction into 2D pixel indices, and then checks whether the pixel
+ * indices fall in the skymap.
+ ***************************************************************************/
+bool GSkymap::contains(const GSkyDir& dir) const
+{
+    // Convert sky direction into sky pixel
+    GSkyPixel pixel = dir2pix(dir);
+    
+    // Return location flag
+    return (contains(pixel));
+}
+
+
+/***********************************************************************//**
+ * @brief Checks if sky map pixel falls in map
+ *
+ * @param[in] pixel Sky map pixel.
+ * @return Trus if pixels is within map, false otherwise.
+ *
+ * Checks whether the specified sky map @p pixel falls within the skymap
+ * or not. 
+ ***************************************************************************/
+bool GSkymap::contains(const GSkyPixel& pixel) const
+{
+    // Initialise containment flag
+    bool inmap = false;
+
+    // Test 1D pixels
+    if (pixel.is1D()) {
+        if (pixel.index()+0.5 >= 0.0 && pixel.index()-0.5 < m_num_pixels) {
+            inmap = true;
+        }
+    }
+
+    // Test 2D pixels
+    else if (pixel.is2D()) {
+
+        // If pixel is in range then set containment flag to true
+        if ((pixel.x()+0.5 >= 0.0 && pixel.x()-0.5 < m_num_x) &&
+            (pixel.y()+0.5 >= 0.0 && pixel.y()-0.5 < m_num_y)) {
+            inmap = true;
+        }
+
+    }
+
+    // Return containment flag
+    return inmap;
+}
+
 /***********************************************************************//**
  * @brief Load skymap from FITS file.
@@ -544,6 +920,9 @@ GSkymap* GSkymap::clone(void) const
  * Loads HEALPix and non HEALPix skymaps. First searches for HEALPix map in
  * FITS file by scanning all HDUs for PIXTYPE=HEALPIX. If no HEALPix map has
  * been found then search load first non-empty image.
+ *
+ * @todo Do we have to restrict a HEALPix map to a BinTable and a WCS map
+ * to a Double precision image???
  ***************************************************************************/
 void GSkymap::load(const std::string& filename)
 {
@@ -568,16 +947,12 @@ void GSkymap::load(const std::string&amp; filename)
         // Get pointer to HDU
         hdu = fits.hdu(extno);
-
+        
         // If PIXTYPE keyword equals "HEALPIX" then load map
-        try {
-            if (hdu->string("PIXTYPE") == "HEALPIX") {
-                read_healpix(static_cast<const GFitsBinTable*>(hdu));
-                loaded = true;
-                break;
-            }
-        }
-        catch (GException::fits_key_not_found &e) {
+        if (hdu->hascard("PIXTYPE") && hdu->string("PIXTYPE") == "HEALPIX") {
+            read_healpix(static_cast<const GFitsTable*>(hdu));
+            loaded = true;
+            break;
         }
     } // endfor: looped over HDUs
@@ -597,7 +972,7 @@ void GSkymap::load(const std::string&amp; filename)
             }
             // Load WCS map
-            read_wcs(static_cast<const GFitsImageDouble*>(hdu));
+            read_wcs(static_cast<const GFitsImage*>(hdu));
             loaded = true;
             break;
@@ -674,14 +1049,10 @@ void GSkymap::read(const GFitsHDU* hdu)
         // Initialize load flag
         bool loaded = false;
-        // Try load as HEALPix map
-        try {
-            if (hdu->string("PIXTYPE") == "HEALPIX") {
-                read_healpix(static_cast<const GFitsBinTable*>(hdu));
-                loaded = true;
-            }
-        }
-        catch (GException::fits_key_not_found &e) {
+        // If PIXTYPE keyword equals "HEALPIX" then load map
+        if (hdu->hascard("PIXTYPE") && hdu->string("PIXTYPE") == "HEALPIX") {
+            read_healpix(static_cast<const GFitsTable*>(hdu));
+            loaded = true;
         }
         // ... otherwise try loading as non HEALPix map
@@ -689,7 +1060,7 @@ void GSkymap::read(const GFitsHDU* hdu)
             // Load only if HDU contains an image
             if (hdu->exttype() == 0) {
-                read_wcs(static_cast<const GFitsImageDouble*>(hdu));
+                read_wcs(static_cast<const GFitsImage*>(hdu));
                 loaded = true;
             }
@@ -741,332 +1112,6 @@ void GSkymap::write(GFits* file) const
 /***********************************************************************//**
- * @brief Returns sky direction of pixel
- *
- * @param[in] pix Pixel number (0,1,...,m_num_pixels).
- *
- * @exception GException::wcs
- *            No valid WCS found.
- *
- * Returns sky direction for a given sky map pixel. This methods works for
- * sky maps that have a 1D (e.g. HEALPix) or 2D pixel indexation scheme.
- * Automatic index conversion is provided so that 1D schemes need only to
- * implement 1D conversion methods while 2D schemes need only to implement
- * 2D conversion methods.
- ***************************************************************************/
-GSkyDir GSkymap::pix2dir(const int& pix) const
-{
-    // Throw error if sky projection is not valid
-    if (m_proj == NULL) {
-        throw GException::wcs(G_PIX2DIR, "No valid sky projection found.");
-    }
-
-    // Determine sky direction from pixel. Use 2D version if sky map is
-    // 2D, otherwise use 1D version.
-    GSkyDir dir = (m_num_x == 0) ? m_proj->pix2dir(pix)
-                                 : m_proj->xy2dir(pix2xy(pix));
-
-    // Return sky direction
-    return dir;
-}
-
-
-/***********************************************************************//**
- * @brief Returns pixel index for a given sky direction
- *
- * @param[in] dir Sky direction.
- *
- * @exception GException::wcs
- *            No valid WCS found.
- *
- * Returns sky map pixel for a given sky direction. This methods works for
- * sky maps that have a 1D (e.g. HEALPix) or 2D pixel indexation scheme.
- * Automatic index conversion is provided so that 1D schemes need only to
- * implement 1D conversion methods while 2D schemes need only to implement
- * 2D conversion methods.
- ***************************************************************************/
-int GSkymap::dir2pix(const GSkyDir& dir) const
-{
-    // Throw error if WCS is not valid
-    if (m_proj == NULL) {
-        throw GException::wcs(G_DIR2PIX, "No valid sky projection found.");
-    }
-
-    // Determine 1D pixel index for a given sky direction
-    int pix = (m_num_x == 0) ? m_proj->dir2pix(dir)
-                             : xy2pix(m_proj->dir2xy(dir));
-
-    // Return pixel index
-    return pix;
-}
-
-
-/***********************************************************************//**
- * @brief Returns sky direction of pixel
- *
- * @param[in] pix Sky map pixel.
- *
- * @exception GException::wcs
- *            No valid WCS found.
- *
- * Returns sky direction for a given sky map pixel. This methods works for
- * sky maps that have a 1D (e.g. HEALPix) or 2D pixel indexation scheme.
- * Automatic index conversion is provided so that 1D schemes need only to
- * implement 1D conversion methods while 2D schemes need only to implement
- * 2D conversion methods.
- ***************************************************************************/
-GSkyDir GSkymap::xy2dir(const GSkyPixel& pix) const
-{
-    // Throw error if WCS is not valid
-    if (m_proj == NULL) {
-        throw GException::wcs(G_XY2DIR, "No valid sky projection found.");
-    }
-
-    // Determine sky direction from pixel. Use 2D version if sky map is
-    // 2D, otherwise use 1D version.
-    GSkyDir dir = (m_num_x == 0) ? m_proj->pix2dir(xy2pix(pix))
-                                 : m_proj->xy2dir(pix);
-
-    // Return sky direction
-    return dir;
-}
-
-
-/***********************************************************************//**
- * @brief Returns sky map pixel for a given sky direction
- *
- * @param[in] dir Sky direction.
- *
- * @exception GException::wcs
- *            No valid WCS found.
- *
- * Returns sky map pixel for a given sky direction. This methods works for
- * sky maps that have a 1D (e.g. HEALPix) or 2D pixel indexation scheme.
- * Automatic index conversion is provided so that 1D schemes need only to
- * implement 1D conversion methods while 2D schemes need only to implement
- * 2D conversion methods.
- ***************************************************************************/
-GSkyPixel GSkymap::dir2xy(const GSkyDir& dir) const
-{
-    // Throw error if WCS is not valid
-    if (m_proj == NULL) {
-        throw GException::wcs(G_DIR2XY, "No valid sky projection found.");
-    }
-
-    // Determine 1D pixel index for a given sky direction
-    GSkyPixel pix = (m_num_x == 0) ? pix2xy(m_proj->dir2pix(dir))
-                                   : m_proj->dir2xy(dir);
-
-    // Return pixel index
-    return pix;
-}
-
-
-/***********************************************************************//**
- * @brief Returns solid angle of pixel
- *
- * @param[in] pix Pixel number (0,1,...,m_num_pixels).
- *
- * @exception GException::wcs
- *            No valid WCS found.
- ***************************************************************************/
-double GSkymap::omega(const int& pix) const
-{
-    // Throw error if WCS is not valid
-    if (m_proj == NULL) {
-        throw GException::wcs(G_OMEGA1, "No valid sky projection found.");
-    }
-
-    // Determine solid angle from pixel. Use 2D version if sky map is
-    // 2D, otherwise use 1D version.
-    double omega = (m_num_x == 0) ? m_proj->omega(pix)
-                                  : m_proj->omega(pix2xy(pix));
-
-    // Return solid angle
-    return omega;
-}
-
-
-/***********************************************************************//**
- * @brief Returns solid angle of pixel
- *
- * @param[in] pix Sky map pixel.
- *
- * @exception GException::wcs
- *            No valid WCS found.
- ***************************************************************************/
-double GSkymap::omega(const GSkyPixel& pix) const
-{
-    // Throw error if WCS is not valid
-    if (m_proj == NULL) {
-        throw GException::wcs(G_OMEGA2, "No valid sky projection found.");
-    }
-
-    // Determine solid angle from pixel. Use 2D version if sky map is
-    // 2D, otherwise use 1D version.
-    double omega = (m_num_x == 0) ? m_proj->omega(xy2pix(pix))
-                                  : m_proj->omega(pix);
-
-    // Return solid angle
-    return omega;
-}
-
-
-/***********************************************************************//**
- * @brief Returns number of pixels
- ***************************************************************************/
-int GSkymap::npix(void) const
-{
-    // Return number of pixels
-    return m_num_pixels;
-}
-
-
-/***********************************************************************//**
- * @brief Returns number of pixels in x coordinate
- ***************************************************************************/
-int GSkymap::nx(void) const
-{
-    // Return number of pixels
-    return m_num_x;
-}
-
-
-/***********************************************************************//**
- * @brief Returns number of pixels in y coordinate
- ***************************************************************************/
-int GSkymap::ny(void) const
-{
-    // Return number of pixels
-    return m_num_y;
-}
-
-
-/***********************************************************************//**
- * @brief Returns number of maps
- ***************************************************************************/
-int GSkymap::nmaps(void) const
-{
-    // Return number of maps
-    return m_num_maps;
-}
-
-
-/***********************************************************************//**
- * @brief Converts 2D index (x,y) into 1D pixel index
- *
- * @param[in] pix 2D pixel index.
- *
- * The (x,y) value is rounded to nearest integers before conversion. The x 
- * axis is assumed as the most rapidely varying index.
- ***************************************************************************/
-int GSkymap::xy2pix(const GSkyPixel& pix) const
-{
-    // Get x and y indices by rounding the (x,y) values
-    int ix = int(pix.x()+0.5);
-    int iy = int(pix.y()+0.5);
-
-    // Return index
-    return (ix+iy*m_num_x);
-}
-
-
-/***********************************************************************//**
- * @brief Converts 1D pixel index into 2D index (x,y)
- *
- * @param[in] pix 1D pixel index.
- *
- * If skymap is not 2D the pixel index is returned in the x element of the
- * GSkyPixel object.
- ***************************************************************************/
-GSkyPixel GSkymap::pix2xy(const int& pix) const
-{
-    // Get x and y indices
-    double x;
-    double y;
-    if (m_num_x != 0) {
-        x = double(pix % m_num_x);
-        y = double(pix / m_num_x);
-    }
-    else {
-        x = double(pix);
-        y = 0.0;
-    }
-
-    // Set pixel
-    GSkyPixel pixel(x, y);
-
-    // Return pixel
-    return pixel;
-}
-
-
-/***********************************************************************//**
- * @brief Set sky projection
- *
- * @param[in] proj Sky projection.
- *
- * Sets the projection from celestial to pixel coordinates. The method
- * performs a deep copy of @p proj, allowing to destroy the argument after
- * using the method. 
- ***************************************************************************/
-void GSkymap::projection(const GSkyProjection& proj)
-{
-    // Free any existing WCS
-    if (m_proj != NULL) delete m_proj;
-
-    // Clone input WCS
-    m_proj = proj.clone();
-
-    // Return
-    return;
-}
-
-
-/***********************************************************************//**
- * @brief Verifies if sky direction falls in map
- *
- * @param[in] dir Sky direction.
- *
- * This method checks if the specified sky direction falls within the pixels
- * covered by the skymap. The method uses the dir2xy method to convert the
- * sky direction into 2D pixel indices, and then checks whether the pixel
- * indices fall in the skymap.
- ***************************************************************************/
-bool GSkymap::contains(const GSkyDir& dir) const
-{
-    // Convert sky direction into sky pixel
-    GSkyPixel pixel = dir2xy(dir);
-    
-    // Return location flag
-    return (contains(pixel));
-}
-
-
-/***********************************************************************//**
- * @brief Verifies if sky pixel falls in map
- *
- * @param[in] pixel Sky pixel.
- *
- * This method checks if the specified sky pixel is within the skymap.
- ***************************************************************************/
-bool GSkymap::contains(const GSkyPixel& pixel) const
-{
-    // Initialise location flag
-    bool inmap = false;
-
-    // If pixel is in range then set location flag to true
-    if ((pixel.x()+0.5 >= 0.0 && pixel.x()-0.5 < m_num_x) &&
-        (pixel.y()+0.5 >= 0.0 && pixel.y()-0.5 < m_num_y)) {
-        inmap = true;
-    }
-
-    // Return location flag
-    return inmap;
-}
-
-
-/***********************************************************************//**
  * @brief Print sky map
  *
  * @param[in] chatter Chattiness (defaults to NORMAL).
@@ -1088,7 +1133,7 @@ std::string GSkymap::print(const GChatter&amp; chatter) const
         result.append(gammalib::str(m_num_pixels));
         result.append("\n"+gammalib::parformat("Number of maps"));
         result.append(gammalib::str(m_num_maps));
-        if (m_proj != NULL && m_proj->code() != "HPX") {
+        if (m_proj != NULL && m_proj->size() == 2) {
             result.append("\n"+gammalib::parformat("X axis dimension"));
             result.append(gammalib::str(m_num_x));
             result.append("\n"+gammalib::parformat("Y axis dimension"));
@@ -1181,7 +1226,7 @@ void GSkymap::copy_members(const GSkymap&amp; map)
     // Copy pixels
     if (size > 0 && map.m_pixels != NULL) {
         alloc_pixels();
-        for (int i = 0; i <  size; ++i) {
+        for (int i = 0; i < size; ++i) {
             m_pixels[i] = map.m_pixels[i];
         }
     }
@@ -436,35 +436,9 @@ void GWcs::write(GFitsHDU* hdu) const
 /***********************************************************************//**
  * @brief Returns solid angle of pixel in units of steradians
  *
- * @param[in] pix Pixel index.
+ * @param[in] pixel Pixel index (x,y)
  *
- * This dummy method throws an error when called as pixel indices are not
- * implemented for a HealPix grid. Maybe we should drop this method and
- * implement a general GSkyPixel class tha also handles HealPix?
- *
- * @todo Think about implementation of sky pixel for HealPix data. We may
- *       use only the first argument, and even carry a usage flag in
- *       GSkyPixel
- ***************************************************************************/
-double GWcs::omega(const int& pix) const
-{
-    // Set error message
-    std::string message = "Method not defined for general WCS projection.";
-
-    // Throw error
-    throw GException::wcs(G_OMEGA, message);
-
-    // Return
-    return 0.0;
-}
-
-
-/***********************************************************************//**
- * @brief Returns solid angle of pixel in units of steradians
- *
- * @param[in] pix Pixel index (x,y)
- *
- * Estimate solid angles of pixel by compuing the coordinates in the 4 pixel
+ * Estimate solid angles of pixel by computing the coordinates in the 4 pixel
  * corners. The surface is computed using a cartesian approximation:
  *           a
  *     1-----5-----2                 a+b
@@ -475,15 +449,15 @@ double GWcs::omega(const int&amp; pix) const
  * rotated sky maps. Something more intelligent should be implemented in
  * the future.
  ***************************************************************************/
-double GWcs::omega(const GSkyPixel& pix) const
+double GWcs::omega(const GSkyPixel& pixel) const
 {
     // Get the sky directions of the 6 points
-    GSkyDir dir1 = xy2dir(GSkyPixel(pix.x()-0.5, pix.y()-0.5));
-    GSkyDir dir2 = xy2dir(GSkyPixel(pix.x()+0.5, pix.y()-0.5));
-    GSkyDir dir3 = xy2dir(GSkyPixel(pix.x()+0.5, pix.y()+0.5));
-    GSkyDir dir4 = xy2dir(GSkyPixel(pix.x()-0.5, pix.y()+0.5));
-    GSkyDir dir5 = xy2dir(GSkyPixel(pix.x(), pix.y()-0.5));
-    GSkyDir dir6 = xy2dir(GSkyPixel(pix.x(), pix.y()+0.5));
+    GSkyDir dir1 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()-0.5));
+    GSkyDir dir2 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()-0.5));
+    GSkyDir dir3 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()+0.5));
+    GSkyDir dir4 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+0.5));
+    GSkyDir dir5 = pix2dir(GSkyPixel(pixel.x(),     pixel.y()-0.5));
+    GSkyDir dir6 = pix2dir(GSkyPixel(pixel.x(),     pixel.y()+0.5));
     // Compute distances between sky directions
     double a = dir1.dist(dir2);
@@ -499,64 +473,16 @@ double GWcs::omega(const GSkyPixel&amp; pix) const
 /***********************************************************************//**
- * @brief Returns sky direction of pixel (DUMMY METHOD)
- *
- * @param[in] pix Pixel number (0,1,...,m_num_pixels).
- *
- * @exception GException::wcs
- *            Method not defined for projection.
- *
- * This method should be overloaded by the appropriate method of the derived
- * class. It thus should never get called.
- ***************************************************************************/
-GSkyDir GWcs::pix2dir(const int& pix) const
-{
-    // Set error message
-    std::string message = "Method not defined for general WCS projection.";
-
-    // Throw error
-    throw GException::wcs(G_PIX2DIR, message);
-
-    // Define direction
-    GSkyDir dir;
-
-    // Return
-    return dir;
-}
-
-
-/***********************************************************************//**
- * @brief Returns pixel for sky direction (DUMMY METHOD)
- *
- * @param[in] dir Sky direction.
- *
- * @exception GException::wcs
- *            Method not defined for projection.
- *
- * This method should be overloaded by the appropriate method of the derived
- * class. It thus should never get called.
- ***************************************************************************/
-int GWcs::dir2pix(const GSkyDir& dir) const
-{
-    // Set error message
-    std::string message = "Method not defined for general WCS projection.";
-
-    // Throw error
-    throw GException::wcs(G_DIR2PIX, message);
-
-    // Return
-    return 0;
-}
-
-
-/***********************************************************************//**
- * @brief Returns sky direction of pixel
+ * @brief Returns sky direction of sky map pixel
  *
- * @param[in] pix Sky pixel.
+ * @param[in] pixel Sky map pixel.
+ * @return Sky direction.
  *
- * Note that pixel indices start from 0.
+ * Returns the sky direction of a sky map @p pixel. Note that the sky
+ * map pixel values start from 0 while the WCS pixel reference starts from
+ * 1.
  ***************************************************************************/
-GSkyDir GWcs::xy2dir(const GSkyPixel& pix) const
+GSkyDir GWcs::pix2dir(const GSkyPixel& pixel) const
 {
     // Allocate memory for transformation
     double pixcrd[2];
@@ -568,8 +494,8 @@ GSkyDir GWcs::xy2dir(const GSkyPixel&amp; pix) const
     // Set sky pixel. We have to add 1.0 here as the WCS pixel reference
     // (CRPIX) starts from one while GSkyPixel starts from 0.
-    pixcrd[0] = pix.x() + 1.0;
-    pixcrd[1] = pix.y() + 1.0;
+    pixcrd[0] = pixel.x() + 1.0;
+    pixcrd[1] = pixel.y() + 1.0;
     // Trasform pixel-to-world coordinate
     wcs_p2s(1, 2, pixcrd, imgcrd, &phi, &theta, world, &stat);
@@ -585,7 +511,7 @@ GSkyDir GWcs::xy2dir(const GSkyPixel&amp; pix) const
     // Debug: Dump transformation steps
     #if defined(G_XY2DIR_DEBUG)
-    std::cout << "xy2dir: pixel=" << pix
+    std::cout << "xy2dir: pixel=" << pixel
               << " (x,y)=(" << pixcrd[0] << "," << pixcrd[1] << ")"
               << " (phi,theta)=(" << phi << "," << theta << ")"
               << " (lng,lat)=(" << world[0] << "," << world[1] << ")" << std::endl;
@@ -597,14 +523,16 @@ GSkyDir GWcs::xy2dir(const GSkyPixel&amp; pix) const
 /***********************************************************************//**
- * @brief Returns pixel of sky direction
+ * @brief Returns sky map pixel of sky direction
  *
  * @param[in] dir Sky direction.
+ * @return Sky map pixel.
  *
- * Note that GSkyPixel indices start from 0 while the WCS pixel reference
- * starts from 1.
+ * Returns the sky map pixel for a given sky direction. Note that the sky
+ * map pixel values start from 0 while the WCS pixel reference starts from
+ * 1.
  ***************************************************************************/
-GSkyPixel GWcs::dir2xy(const GSkyDir& dir) const
+GSkyPixel GWcs::dir2pix(const GSkyDir& dir) const
 {
     // Allocate memory for transformation
     double pixcrd[2];
@@ -106,10 +106,10 @@ double TestGSky::wcs_forth_back_pixel(GWcs* wcs, int nx, int ny, double&amp; crpix1,
             GSkyPixel pix_in(x,y);
             // Forth: Transform to world
-            GSkyDir dir = wcs->xy2dir(pix_in);
+            GSkyDir dir = wcs->pix2dir(pix_in);
             // Back: Transform to pixel
-            GSkyPixel pix_out = wcs->dir2xy(dir);
+            GSkyPixel pix_out = wcs->dir2pix(dir);
             // Compute distance
             double dx      = pix_in.x()-pix_out.x();
@@ -177,7 +177,7 @@ double TestGSky::wcs_copy(GWcs* wcs, int nx, int ny, double&amp; crpix1, double&amp; crp
         for (int iy = -ny; iy <= ny; ++iy) {
             // Set y value
-            double    y = double(iy) + crpix2;
+            double y = double(iy) + crpix2;
             // Skip pixels outside valid range (they are not expected to
             // transform bijectively)
@@ -189,8 +189,8 @@ double TestGSky::wcs_copy(GWcs* wcs, int nx, int ny, double&amp; crpix1, double&amp; crp
             GSkyPixel pix_in(x,y);
             // Transform to world
-            GSkyDir dir1 = wcs->xy2dir(pix_in);
-            GSkyDir dir2 = cpy->xy2dir(pix_in);
+            GSkyDir dir1 = wcs->pix2dir(pix_in);
+            GSkyDir dir2 = cpy->pix2dir(pix_in);
             double angle = dir1.dist_deg(dir2);
             if (angle > angle_max) {
                 angle_max = angle;
@@ -207,8 +207,8 @@ double TestGSky::wcs_copy(GWcs* wcs, int nx, int ny, double&amp; crpix1, double&amp; crp
             #endif
             // Transform to pixel
-            GSkyPixel pix_out1 = wcs->dir2xy(dir1);
-            GSkyPixel pix_out2 = cpy->dir2xy(dir2);
+            GSkyPixel pix_out1 = wcs->dir2pix(dir1);
+            GSkyPixel pix_out2 = cpy->dir2pix(dir2);
             // Compute distance
             double dx   = pix_out1.x()-pix_out2.x();
@@ -715,8 +715,8 @@ void TestGSky::test_GSkymap_wcs_construct(void)
         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);
+                GSkyPixel pixel    = map1.dir2pix(dir);
+                GSkyDir   dir_back = map1.pix2dir(pixel);
                 double    dist     = dir.dist_deg(dir_back);
                 if (dist > eps) {
                     throw exception_failure("Sky direction differs: dir="+dir.print()+" pixel="+pixel.print()+" dir_back"+ dir_back.print()+" dist="+gammalib::str(dist)+" deg");
@@ -73,10 +73,11 @@ class Test(GPythonTestSuite):
         # Append tests
         self.append(self.test_skymap_healpix, "Test HEALPix map")
         self.append(self.test_skymap_ait, "Test AIT projection map")
+        self.append(self.test_skymap_azp, "Test AZP projection map")
         self.append(self.test_skymap_car, "Test CAR projection map")
-        self.append(self.test_skymap_tan, "Test TAN projection map")
+        self.append(self.test_skymap_mer, "Test MER projection map")
         self.append(self.test_skymap_stg, "Test STG projection map")
-        self.append(self.test_skymap_azp, "Test AZP projection map")
+        self.append(self.test_skymap_tan, "Test TAN projection map")
         self.append(self.test_fk5_to_galactic, "Test FK5 to Galactic coordinate conversion")
         # Return
@@ -87,21 +88,21 @@ class Test(GPythonTestSuite):
         """
         Control coordinate and pixel transformation methods.
         """
-        # Control pix2dir / dir2pix methods for all pixels
+        # Control inx2dir/dir2inx methods for all pixels
         for i in range(pixels.npix()):
-            dir = pixels.pix2dir(i)
-            ipix = pixels.dir2pix(dir)
-            msg = map + " pix2dir/dir2pix check for pixel " + str(i)
-            err = map + " GSkymap trouble with pixel " + str(i) + " (" + str(ipix) + \
+            dir = pixels.inx2dir(i)
+            inx = pixels.dir2inx(dir)
+            msg = map + " inx2dir/dir2inx check for pixel " + str(i)
+            err = map + " GSkymap trouble with pixel " + str(i) + " (" + str(inx) + \
                 "), RA=" + str(dir.ra() * 180 / pi) + ", Dec=" + str(dir.dec() * 180 / pi)
-            self.test_assert(i == ipix, msg, err)
+            self.test_assert(i == inx, msg, err)
         # Control SkyDir coordinate transformation for all pixels
         for i in range(pixels.npix()):
-            dir = pixels.pix2dir(i)
+            dir     = pixels.inx2dir(i)
             dir_new = GSkyDir()
             dir_new.lb(dir.l(), dir.b())
-            dra = abs(dir.ra() - dir_new.ra())
+            dra     = abs(dir.ra() - dir_new.ra())
             if (dra >= 5.0):
                 dra -= 2.0 * pi
             ddec = abs(dir.dec() - dir_new.dec())
@@ -196,6 +197,17 @@ class Test(GPythonTestSuite):
         # Return
         return
+    # Test AZP projection
+    def test_skymap_azp(self):
+        """
+        Test AZP projection.
+        """
+        # Execute generic test
+        self.test_skymap_proj("AZP")
+
+        # Return
+        return
+
     # Test CAR projection
     def test_skymap_car(self):
         """
@@ -207,13 +219,13 @@ class Test(GPythonTestSuite):
         # Return
         return
-    # Test TAN projection
-    def test_skymap_tan(self):
+    # Test MER projection
+    def test_skymap_mer(self):
         """
-        Test TAN projection.
+        Test MER projection.
         """
         # Execute generic test
-        self.test_skymap_proj("TAN")
+        self.test_skymap_proj("MER")
         # Return
         return
@@ -229,13 +241,13 @@ class Test(GPythonTestSuite):
         # Return
         return
-    # Test AZP projection
-    def test_skymap_azp(self):
+    # Test TAN projection
+    def test_skymap_tan(self):
         """
-        Test AZP projection.
+        Test TAN projection.
         """
         # Execute generic test
-        self.test_skymap_proj("AZP")
+        self.test_skymap_proj("TAN")
         # Return
         return