diff --git a/include/GSkymap.hpp b/include/GSkymap.hpp index fe2931e..9a52db0 100644 --- a/include/GSkymap.hpp +++ b/include/GSkymap.hpp @@ -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) diff --git a/include/GWcs.hpp b/include/GWcs.hpp index 61fc99d..ab55385 100644 --- a/include/GWcs.hpp +++ b/include/GWcs.hpp @@ -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 diff --git a/include/GWcsCAR.hpp b/include/GWcsCAR.hpp index e3a7a77..76d5108 100644 --- a/include/GWcsCAR.hpp +++ b/include/GWcsCAR.hpp @@ -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; diff --git a/include/GWcsHPX.hpp b/include/GWcsHPX.hpp index 6f59646..7d3fac1 100644 --- a/include/GWcsHPX.hpp +++ b/include/GWcsHPX.hpp @@ -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; diff --git a/src/sky/GSkymap.cpp b/src/sky/GSkymap.cpp index 0ff8391..3087fa1 100644 --- a/src/sky/GSkymap.cpp +++ b/src/sky/GSkymap.cpp @@ -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& wcs, const std::string& 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& 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& 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& 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 = diff --git a/src/sky/GWcs.cpp b/src/sky/GWcs.cpp index 8e89ed6..2d0965c 100644 --- a/src/sky/GWcs.cpp +++ b/src/sky/GWcs.cpp @@ -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& 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& GWcs::operator= (const GWcs& 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& 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; diff --git a/src/sky/GWcsCAR.cpp b/src/sky/GWcsCAR.cpp index f1d5005..22ee48b 100644 --- a/src/sky/GWcsCAR.cpp +++ b/src/sky/GWcsCAR.cpp @@ -188,12 +188,46 @@ GWcsCAR& GWcsCAR::operator= (const GWcsCAR& 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& pix) const void GWcsCAR::init_members(void) { // Initialise members + m_type = "CAR"; // Return return; @@ -315,7 +350,7 @@ std::ostream& operator<< (std::ostream& os, const GWcsCAR& wcs) os << "=== GWcsCAR ===" << std::endl; // Add WCS information - wcs.dump_wcs(os); + wcs.wcs_dump(os); // Return output stream return os; diff --git a/src/sky/GWcsHPX.cpp b/src/sky/GWcsHPX.cpp index 25f698f..f27b676 100644 --- a/src/sky/GWcsHPX.cpp +++ b/src/sky/GWcsHPX.cpp @@ -211,16 +211,6 @@ GWcsHPX& GWcsHPX::operator= (const GWcsHPX& 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& ordering) void GWcsHPX::init_members(void) { // Initialise members + m_type = "HPX"; m_nside = 0; m_npface = 0; m_ncap = 0; diff --git a/test/test_GSky.cpp b/test/test_GSky.cpp index b2e9064..c6e0cd3 100644 --- a/test/test_GSky.cpp +++ b/test/test_GSky.cpp @@ -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 << map1 << 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;