lsst::afw::image::Wcs Class Reference

Implementation of the WCS standard for a any projection. More...

#include <Wcs.h>

Inheritance diagram for lsst::afw::image::Wcs:

Public Types
typedef boost::shared_ptr< Wcs >	Ptr
typedef boost::shared_ptr< Wcs const >	ConstPtr
Public Types inherited from lsst::daf::base::Persistable
typedef boost::shared_ptr < Persistable >	Ptr
Public Types inherited from lsst::daf::base::Citizen
enum	{ magicSentinel = 0xdeadbeef }
typedef unsigned long	memId
	Type of the block's ID. More...
typedef memId(*	memNewCallback )(const memId cid)
	A function used to register a callback. More...
typedef memId(*	memCallback )(const Citizen *ptr)

Public Member Functions
	Wcs (lsst::afw::geom::Point2D const &crval, lsst::afw::geom::Point2D const &crpix, Eigen::Matrix2d const &CD, std::string const &ctype1="RA---TAN", std::string const &ctype2="DEC--TAN", double equinox=2000, std::string const &raDecSys="ICRS", std::string const &cunits1="deg", std::string const &cunits2="deg")
	Create a Wcs object with some known information. More...
virtual	~Wcs ()
virtual Ptr	clone (void) const
bool	operator== (Wcs const &other) const
bool	operator!= (Wcs const &other) const
lsst::afw::coord::Coord::Ptr	getSkyOrigin () const
	Returns CRVAL. More...
lsst::afw::geom::Point2D	getPixelOrigin () const
	Returns CRPIX (corrected to LSST convention). More...
Eigen::Matrix2d	getCDMatrix () const
	Returns CD matrix. You would never have guessed that from the name. More...
virtual void	flipImage (int flipLR, int flipTB, lsst::afw::geom::Extent2I dimensions) const
	Flip CD matrix around the y-axis. More...
virtual void	rotateImageBy90 (int nQuarter, lsst::afw::geom::Extent2I dimensions) const
virtual boost::shared_ptr < lsst::daf::base::PropertyList >	getFitsMetadata () const
	Return a PropertyList containing FITS header keywords that can be used to save the Wcs.x. More...
bool	isFlipped () const
double	pixArea (lsst::afw::geom::Point2D pix00) const
	Sky area covered by a pixel at position pix00 in units of square degrees. More...
geom::Angle	pixelScale () const
	Returns the pixel scale [Angle/pixel]. More...
boost::shared_ptr< coord::Coord >	pixelToSky (double pix1, double pix2) const
	Convert from celestial coordinates to pixel coordinates. More...
boost::shared_ptr< coord::Coord >	pixelToSky (lsst::afw::geom::Point2D const &pixel) const
	Convert from celestial coordiantes to pixel coordinates. More...
void	pixelToSky (double pixel1, double pixel2, geom::Angle &sky1, geom::Angle &sky2) const
	Convert from celestial coordiantes to pixel coordinates. More...
geom::Point2D	skyToPixel (geom::Angle sky1, geom::Angle sky2) const
	Convert from sky coordinates (e.g ra/dec) to pixel positions. More...
geom::Point2D	skyToPixel (coord::Coord const &coord) const
	Convert from sky coordinates (e.g ra/dec) to pixel positions. More...
geom::Point2D	skyToIntermediateWorldCoord (coord::Coord const &coord) const
	Convert from sky coordinates (e.g ra/dec) to intermediate world coordinates. More...
virtual bool	hasDistortion () const
geom::LinearTransform	getLinearTransform () const
geom::AffineTransform	linearizePixelToSky (coord::Coord const &coord, geom::AngleUnit skyUnit=geom::degrees) const
	Return the local linear approximation to Wcs::pixelToSky at a point given in sky coordinates. More...
geom::AffineTransform	linearizePixelToSky (geom::Point2D const &pix, geom::AngleUnit skyUnit=geom::degrees) const
	Return the local linear approximation to Wcs::pixelToSky at a point given in pixel coordinates. More...
geom::AffineTransform	linearizeSkyToPixel (coord::Coord const &coord, geom::AngleUnit skyUnit=geom::degrees) const
	Return the local linear approximation to Wcs::skyToPixel at a point given in sky coordinates. More...
geom::AffineTransform	linearizeSkyToPixel (geom::Point2D const &pix, geom::AngleUnit skyUnit=geom::degrees) const
	Return the local linear approximation to Wcs::skyToPixel at a point given in pixel coordinates. More...
void	shiftReferencePixel (geom::Extent2D const &d)
void	shiftReferencePixel (double dx, double dy)
	Move the pixel reference position by (dx, dy) Used when persisting and retrieving sub-images. The lsst convention is that Wcs returns pixel position (which is based on position in the parent image), but the fits convention is to return pixel index (which is bases on position in the sub-image). In order that the fits files we create make sense to other fits viewers, we change to the fits convention when writing out images. More...
virtual bool	isPersistable () const
	Whether the Wcs is persistable using afw::table::io archives. More...
Public Member Functions inherited from lsst::daf::base::Persistable
	Persistable (void)
virtual	~Persistable (void)
template<class Archive >
void	serialize (Archive &, unsigned int const)
Public Member Functions inherited from lsst::daf::base::Citizen
	Citizen (const std::type_info &)
	Citizen (Citizen const &)
	~Citizen ()
Citizen &	operator= (Citizen const &)
std::string	repr () const
	Return a string representation of a Citizen. More...
void	markPersistent (void)
	Mark a Citizen as persistent and not destroyed until process end. More...
memId	getId () const
	Return the Citizen's ID. More...
Public Member Functions inherited from lsst::afw::table::io::Persistable
void	writeFits (std::string const &fileName, std::string const &mode="w") const
	Write the object to a regular FITS file. More...
void	writeFits (fits::MemFileManager &manager, std::string const &mode="w") const
	Write the object to a FITS image in memory. More...
void	writeFits (fits::Fits &fitsfile) const
	Write the object to an already-open FITS object. More...
virtual	~Persistable ()

Protected Member Functions
virtual std::string	getPersistenceName () const
	Return the unique name used to persist this object and look up its factory. More...
virtual std::string	getPythonModule () const
	Return the fully-qualified Python module that should be imported to guarantee that its factory is registered. More...
virtual void	write (OutputArchiveHandle &handle) const
	Write the object to one or more catalogs. More...
virtual bool	_isSubset (Wcs const &other) const
	Wcs ()
	Construct an invalid Wcs given no arguments. More...
	Wcs (boost::shared_ptr< lsst::daf::base::PropertySet const > const &fitsMetadata)
	Wcs (afw::table::BaseRecord const &record)
	Wcs (Wcs const &rhs)
	Copy constructor. More...
Wcs &	operator= (const Wcs &)
afw::coord::Coord::Ptr	makeCorrectCoord (geom::Angle sky0, geom::Angle sky1) const
	Given a sky position, use the values stored in ctype and radesys to return the correct sub-class of Coord. More...
afw::coord::Coord::Ptr	convertCoordToSky (coord::Coord const &coord) const
virtual geom::AffineTransform	linearizePixelToSkyInternal (geom::Point2D const &pix, coord::Coord const &coord, geom::AngleUnit skyUnit) const
virtual geom::AffineTransform	linearizeSkyToPixelInternal (geom::Point2D const &pix, coord::Coord const &coord, geom::AngleUnit skyUnit) const
void	initWcsLibFromFits (boost::shared_ptr< lsst::daf::base::PropertySet const > const &fitsMetadata)
	Parse a fits header, extract the relevant metadata and create a Wcs object. More...
void	_initWcs ()
void	_setWcslibParams ()
Protected Member Functions inherited from lsst::afw::table::io::Persistable
	Persistable ()
	Persistable (Persistable const &other)
void	operator= (Persistable const &other)

Protected Attributes
struct wcsprm *	_wcsInfo
int	_nWcsInfo
int	_relax
	Degree of permissiveness for wcspih (0 for strict); see wcshdr.h for details. More...
int	_wcsfixCtrl
	Do potentially unsafe translations of non-standard unit strings? 0/1 = no/yes. More...
int	_wcshdrCtrl
	Controls messages to stderr from wcshdr (0 for none); see wcshdr.h for details. More...
int	_nReject
coord::CoordSystem	_coordSystem

Private Member Functions
void	initWcsLib (geom::Point2D const &crval, geom::Point2D const &crpix, Eigen::Matrix2d const &CD, std::string const &ctype1, std::string const &ctype2, double equinox, std::string const &raDecSys, std::string const &cunits1, std::string const &cunits2)
	Manually initialise a wcs struct using values passed by the constructor. More...
virtual void	pixelToSkyImpl (double pixel1, double pixel2, geom::Angle skyTmp[2]) const
virtual geom::Point2D	skyToPixelImpl (geom::Angle sky1, geom::Angle sky2) const

Friends
class	WcsFactory
Wcs::Ptr	makeWcs (boost::shared_ptr< lsst::daf::base::PropertySet > const &fitsMetadata, bool stripMetadata)
	Create a Wcs of the correct class using a fits header. More...

Additional Inherited Members
Static Public Member Functions inherited from lsst::daf::base::Citizen
static bool	hasBeenCorrupted ()
	Check all allocated blocks for corruption. More...
static memId	getNextMemId ()
	Return the memId of the next object to be allocated. More...
static int	init ()
	Called once when the memory system is being initialised. More...
static int	census (int, memId startingMemId=0)
	How many active Citizens are there? More...
static void	census (std::ostream &stream, memId startingMemId=0)
	Print a list of all active Citizens to stream, sorted by ID. More...
static const std::vector < const Citizen * > *	census ()
	Return a (newly allocated) std::vector of active Citizens sorted by ID. More...
static memId	setNewCallbackId (memId id)
	Call the NewCallback when block is allocated. More...
static memId	setDeleteCallbackId (memId id)
	Call the current DeleteCallback when block is deleted. More...
static memNewCallback	setNewCallback (memNewCallback func)
	Set the NewCallback function. More...
static memCallback	setDeleteCallback (memCallback func)
	Set the DeleteCallback function. More...
static memCallback	setCorruptionCallback (memCallback func)
	Set the CorruptionCallback function. More...
Static Public Member Functions inherited from lsst::afw::table::io::PersistableFacade< Wcs >
static boost::shared_ptr< Wcs >	readFits (fits::Fits &fitsfile)
	Read an object from an already open FITS object. More...
static boost::shared_ptr< Wcs >	readFits (std::string const &fileName, int hdu=0)
	Read an object from a regular FITS file. More...
static boost::shared_ptr< Wcs >	readFits (fits::MemFileManager &manager, int hdu=0)
	Read an object from a FITS file in memory. More...
Protected Types inherited from lsst::afw::table::io::Persistable
typedef io::OutputArchiveHandle	OutputArchiveHandle

Detailed Description

Implementation of the WCS standard for a any projection.

Implements a single representation of the World Coordinate System of a two dimensional image The standard is defined in two papers

• Greisen & Calabretta, 2002 A&A 395, 1061
• Calabretta & Greisen, 2002, A&A 395, 1077

In its simplest sense, Wcs is used to convert from position in the sky (in right ascension and declination) to pixel position on an image (and back again). It is, however, much more general than that and can understand a myriad of different coordinate systems.

A wcs can be constructed from a reference position (crval, crpix) and a translation matrix. Alternatively, if you have the header from a fits file, you can create a Wcs object with the makeWcs() function. This function determines whether your Wcs is one the subset of projection systems that is dealt with specially by Lsst, and creates an object of the correct class. Otherwise, a pointer to a Wcs object is returned. Most astronomical images use tangent plane projection, so makeWcs() returns a TanWcs object pointer

*  import lsst.afw.image as afwImg
*  fitsHeader = afwImg.readMetadata(filename)
*  
*  if 0:
*      #This doesn't work
*      wcs = afwImg.Wcs(fitsHeader)
*      
*  wcs = afwImg.makeWcs(fitsHeader)
*  
*  pixelPosition = wcs.skyToPixel(ra, dec)
*  skyPosition = wcs.skyToPixel(xPosition, yPosition)
*  

o[ This class is implemented in by calls to the wcslib library by Mark Calabretta http://www.atnf.csiro.au/people/mcalabre/WCS/

Note that we violate the Wcs standard in one minor way. The standard states that none of the CRPIX or CRVAL keywords are required, for the header to be valid, and the appropriate values should be set to 0.0 if the keywords are absent. This is a recipe for painful bugs in analysis, so we violate the standard by insisting that the keywords CRPIX[1,2] and CRVAL[1,2] are present when reading a header (keywords CRPIX1a etc are also accepted)

Definition at line 107 of file Wcs.h.

Member Typedef Documentation

typedef boost::shared_ptr<Wcs const> lsst::afw::image::Wcs::ConstPtr

Definition at line 114 of file Wcs.h.

typedef boost::shared_ptr<Wcs> lsst::afw::image::Wcs::Ptr

Definition at line 113 of file Wcs.h.

Constructor & Destructor Documentation

Wcs::Wcs	(	lsst::afw::geom::Point2D const &	crval,
		lsst::afw::geom::Point2D const &	crpix,
		Eigen::Matrix2d const &	CD,
		std::string const &	ctype1 = "RA---TAN",
		std::string const &	ctype2 = "DEC--TAN",
		double	equinox = 2000,
		std::string const &	raDecSys = "ICRS",
		std::string const &	cunits1 = "deg",
		std::string const &	cunits2 = "deg"
	)

Create a Wcs object with some known information.

Parameters

crval	The sky position of the reference point
crpix	The pixel position corresponding to crval in LSST units
CD	Matrix describing transformations from pixel to sky positions
ctype1	Projection system used (see description of Wcs)
ctype2	Projection system used (see description of Wcs)
equinox	Equinox of coordinate system, eg 2000 (Julian) or 1950 (Besselian)
raDecSys	System used to describe right ascension or declination, e.g FK4, FK5 or ICRS
cunits1	Units of sky position. One of deg, arcmin or arcsec
cunits2	Units of sky position. One of deg, arcmin or arcsec

Note

LSST units are zero indexed while FITs units are 1 indexed. So a value of crpix stored in a fits header of 127,127 corresponds to a pixel position in LSST units of 128, 128

Definition at line 148 of file Wcs.cc.

 :
    daf::base::Citizen(typeid(this)),
    _wcsInfo(NULL), 
    _nWcsInfo(0), 
    _relax(0), 
    _wcsfixCtrl(0), 
    _wcshdrCtrl(0),
    _nReject(0),
    _coordSystem(static_cast<afwCoord::CoordSystem>(-1))
{
    _setWcslibParams();
    initWcsLib(crval, crpix, CD, 
               ctype1, ctype2,
               equinox, raDecSys,
               cunits1, cunits2);
    _initWcs();
}

Wcs::~Wcs ( )

virtual

Definition at line 555 of file Wcs.cc.

          {
    if (_wcsInfo != NULL) {
        wcsvfree(&_nWcsInfo, &_wcsInfo);
    }
}

lsst::afw::image::Wcs::Wcs ( )

protected

Construct an invalid Wcs given no arguments.

Definition at line 87 of file Wcs.cc.

                     :
    daf::base::Citizen(typeid(this)),
    _wcsInfo(NULL), _nWcsInfo(0), _relax(0), _wcsfixCtrl(0), _wcshdrCtrl(0), _nReject(0),
    _coordSystem(static_cast<afwCoord::CoordSystem>(-1)) {
    _setWcslibParams();
    _initWcs();    
}

Wcs::Wcs ( boost::shared_ptr< lsst::daf::base::PropertySet const > const &  fitsMetadata )

protected

Create a Wcs from a fits header. Don't call this directly. Use makeWcs() instead, which will figure out which (if any) sub-class of Wcs is appropriate

Definition at line 98 of file Wcs.cc.

                                                             :
    daf::base::Citizen(typeid(this)),
    _wcsInfo(NULL), 
    _nWcsInfo(0), 
    _relax(0), 
    _wcsfixCtrl(0), 
    _wcshdrCtrl(0),
    _nReject(0),
    _coordSystem(static_cast<afwCoord::CoordSystem>(-1))
{
    _setWcslibParams();

    initWcsLibFromFits(fitsMetadata);
    _initWcs();
}

lsst::afw::image::Wcs::Wcs ( afw::table::BaseRecord const &  record )

explicitprotected

Definition at line 1153 of file Wcs.cc.

                                          :
    daf::base::Citizen(typeid(this)),
    _wcsInfo(NULL),
    _nWcsInfo(0),
    _relax(0),
    _wcsfixCtrl(0),
    _wcshdrCtrl(0),
    _nReject(0),
    _coordSystem(static_cast<afw::coord::CoordSystem>(-1))
{
    WcsPersistenceHelper const & keys = WcsPersistenceHelper::get();
    if (!record.getSchema().contains(keys.schema)) {
        throw LSST_EXCEPT(
            afw::table::io::MalformedArchiveError,
            "Incorrect schema for Wcs persistence"
        );
    }
    _setWcslibParams();
    Eigen::Matrix2d cd = Eigen::Map<Eigen::Matrix2d const>(record[keys.cd].getData());
    initWcsLib(
        record.get(keys.crval), record.get(keys.crpix), cd,
        record.get(keys.ctype1), record.get(keys.ctype2),
        record.get(keys.equinox), record.get(keys.radesys),
        record.get(keys.cunit1), record.get(keys.cunit2)
    );
    _initWcs();
}

Wcs::Wcs ( afwImg::Wcs const &  rhs )

protected

Copy constructor.

Definition at line 464 of file Wcs.cc.

                              : 
    daf::base::Citizen(typeid(this)),
    _wcsInfo(NULL), 
    _nWcsInfo(rhs._nWcsInfo), 
    _relax(rhs._relax), 
    _wcsfixCtrl(rhs._wcsfixCtrl), 
    _wcshdrCtrl(rhs._wcshdrCtrl),
    _nReject(rhs._nReject),
    _coordSystem(static_cast<afwCoord::CoordSystem>(-1))
{
    
    if (rhs._nWcsInfo > 0) {
        _wcsInfo = static_cast<struct wcsprm *>(calloc(rhs._nWcsInfo, sizeof(struct wcsprm)));
        if (_wcsInfo == NULL) {
            throw LSST_EXCEPT(lsst::pex::exceptions::MemoryError, "Cannot allocate WCS info");
        }

        _wcsInfo->flag = -1;
        int alloc = 1;                  //Unconditionally allocate memory when calling
        for (int i = 0; i != rhs._nWcsInfo; ++i) {
            int status = wcscopy(alloc, &rhs._wcsInfo[i], &_wcsInfo[i]);
            if (status != 0) {
                wcsvfree(&i, &_wcsInfo);
                throw LSST_EXCEPT(lsst::pex::exceptions::MemoryError,
                                  (boost::format("Could not copy WCS: wcscopy status = %d : %s") %
                                   status % wcs_errmsg[status]).str());
            }
        }
    }
    _initWcs();
}

Member Function Documentation

void Wcs::_initWcs ( )

protected

Definition at line 117 of file Wcs.cc.

{
    if (_wcsInfo) {
        _coordSystem = afwCoord::makeCoordEnum(_wcsInfo->radesys);

        // tell WCSlib that values have been updated
        _wcsInfo->flag = 0;
        // and then tell it to do its internal magic.
        int status = wcsset(_wcsInfo);
        if (status != 0) {
            throw LSST_EXCEPT(except::RuntimeError,
                              (boost::format("Failed to setup wcs structure with wcsset. Status %d: %s") %
                               status % wcs_errmsg[status] ).str());
        }
    }
}

bool Wcs::_isSubset ( Wcs const &  other ) const

protectedvirtual

Reimplemented in lsst::afw::image::TanWcs.

Definition at line 530 of file Wcs.cc.

                                         {
    CHECK_NULLS(_wcsInfo, rhs._wcsInfo);
    CHECK_NULLS(_wcsInfo->crval, rhs._wcsInfo->crval);
    CHECK_NULLS(_wcsInfo->cd, rhs._wcsInfo->cd);
    CHECK_NULLS(_wcsInfo->crpix, rhs._wcsInfo->crpix);
    CHECK_NULLS(_wcsInfo->cunit, rhs._wcsInfo->cunit);
    CHECK_NULLS(_wcsInfo->ctype, rhs._wcsInfo->ctype);
    return _nWcsInfo == rhs._nWcsInfo &&
        _coordSystem == rhs._coordSystem &&
        _wcsInfo->naxis == rhs._wcsInfo->naxis &&
        _wcsInfo->equinox == rhs._wcsInfo->equinox &&
        _wcsInfo->altlin == rhs._wcsInfo->altlin &&
        compareArrays(_wcsInfo->crval, rhs._wcsInfo->crval, 2) &&
        compareArrays(_wcsInfo->crpix, rhs._wcsInfo->crpix, 2) &&
        compareArrays(_wcsInfo->cd, rhs._wcsInfo->cd, 4) &&
        compareStringArrays(_wcsInfo->cunit, rhs._wcsInfo->cunit, 2) &&
        compareStringArrays(_wcsInfo->ctype, rhs._wcsInfo->ctype, 2) &&
        skyToPixel(_wcsInfo->crval[0] * afwGeom::degrees,
                   _wcsInfo->crval[1] * afwGeom::degrees) ==
        rhs.skyToPixel(_wcsInfo->crval[0] * afwGeom::degrees,
                       _wcsInfo->crval[1] * afwGeom::degrees) &&
        *pixelToSky(_wcsInfo->crpix[0], _wcsInfo->crpix[1]) ==
        *rhs.pixelToSky(_wcsInfo->crpix[0], _wcsInfo->crpix[1]);
}

void lsst::afw::image::Wcs::_setWcslibParams ( )

protected

Definition at line 68 of file Wcs.cc.

{
    _wcsfixCtrl =                       // ctrl for wcsfix
        2;                              // Translate "H" to "h"
    _wcshdrCtrl =                       // ctrl for wcspih
        2;                              // Report each rejected keyrecord and the reason why it was rejected
    _relax =                            // relax parameter for wcspih;
        WCSHDR_all;                     // Accept all extensions recognized by the parser
}

Wcs::Ptr Wcs::clone ( void  ) const

virtual

Reimplemented in lsst::afw::image::TanWcs.

Definition at line 562 of file Wcs.cc.

                            {
    return Wcs::Ptr(new Wcs(*this));
}

afwCoord::Coord::Ptr Wcs::convertCoordToSky ( coord::Coord const &  coord ) const

protected

Given a Coord (as a shared pointer), return the sky position in the correct coordinate system for this Wcs.

Definition at line 798 of file Wcs.cc.

                                                        {
    return coord.convert(_coordSystem);
}

void Wcs::flipImage ( int  flipLR, int  flipTB, lsst::afw::geom::Extent2I  dimensions  ) const

virtual

Flip CD matrix around the y-axis.

Definition at line 605 of file Wcs.cc.

                                                                            {
    assert(_wcsInfo);
    
    int const naxis = _wcsInfo->naxis;

    //If naxis != 2, I'm not sure if any of what follows is correct
    assert(naxis == 2);
    if (flipLR) {
        _wcsInfo->cd[0] = -_wcsInfo->cd[0];
        _wcsInfo->cd[2] = -_wcsInfo->cd[2];
        _wcsInfo->crpix[0] = -_wcsInfo->crpix[0] + dimensions.getX();
    }
    if (flipTB) {
        _wcsInfo->cd[1] = -_wcsInfo->cd[1];
        _wcsInfo->cd[3] = -_wcsInfo->cd[3];
        _wcsInfo->crpix[1] = -_wcsInfo->crpix[1]+dimensions.getY();
    }

    // tells libwcs to invalidate cached data, since transformation has been modified
    _wcsInfo->flag = 0;
}

Eigen::Matrix2d Wcs::getCDMatrix

(

)

const

Returns CD matrix. You would never have guessed that from the name.

Return the CD matrix.

Definition at line 587 of file Wcs.cc.

                                     {
    assert(_wcsInfo);
    int const naxis = _wcsInfo->naxis;

    //If naxis != 2, I'm not sure if any of what follows is correct
    assert(naxis == 2);
    
    Eigen::Matrix2d C;

    for (int i=0; i< naxis; ++i){
        for (int j=0; j<naxis; ++j) {
            C(i,j) = _wcsInfo->cd[ (i*naxis) + j ];
        }
    }

    return C;
}

PropertyList::Ptr Wcs::getFitsMetadata ( ) const

virtual

Return a PropertyList containing FITS header keywords that can be used to save the Wcs.x.

Return the Wcs as a fits header.

Reimplemented in lsst::afw::image::TanWcs.

Definition at line 678 of file Wcs.cc.

                                           {
    return lsst::afw::formatters::WcsFormatter::generatePropertySet(*this);
}

lsst::afw::geom::LinearTransform Wcs::getLinearTransform

(

)

const

Return the linear part of the Wcs, the CD matrix in FITS speak, as an AffineTransform

See Also

Definition at line 1058 of file Wcs.cc.

{
    return lsst::afw::geom::LinearTransform(getCDMatrix());
}

std::string lsst::afw::image::Wcs::getPersistenceName ( ) const

protectedvirtual

Return the unique name used to persist this object and look up its factory.

Must be less than ArchiveIndexSchema::MAX_NAME_LENGTH characters.

Reimplemented from lsst::afw::table::io::Persistable.

Reimplemented in lsst::afw::image::TanWcs.

Definition at line 1124 of file Wcs.cc.

{ return getWcsPersistenceName(); }

GeomPoint Wcs::getPixelOrigin

(

)

const

Returns CRPIX (corrected to LSST convention).

Return crpix in the lsst convention. Note that this need not be the centre of the image.

Definition at line 577 of file Wcs.cc.

                                    {
    assert(_wcsInfo);
    //Convert from fits units back to lsst units
    double p1 = _wcsInfo->crpix[0] + fitsToLsstPixels;
    double p2 = _wcsInfo->crpix[1] + fitsToLsstPixels;
    return afwGeom::Point2D(p1, p2);
}

std::string lsst::afw::image::Wcs::getPythonModule ( ) const

protectedvirtual

Return the fully-qualified Python module that should be imported to guarantee that its factory is registered.

Must be less than ArchiveIndexSchema::MAX_MODULE_LENGTH characters.

Will be ignored if empty.

Reimplemented from lsst::afw::table::io::Persistable.

Definition at line 1126 of file Wcs.cc.

{ return "lsst.afw.image"; }

CoordPtr Wcs::getSkyOrigin

(

)

const

Returns CRVAL.

Return crval. Note that this need not be the centre of the image.

Definition at line 571 of file Wcs.cc.

                                 {
    assert(_wcsInfo);
    return makeCorrectCoord(_wcsInfo->crval[0] * afwGeom::degrees, _wcsInfo->crval[1] * afwGeom::degrees);
}

virtual bool lsst::afw::image::Wcs::hasDistortion ( ) const

inlinevirtual

Reimplemented in lsst::afw::image::TanWcs.

Definition at line 198 of file Wcs.h.

{    return false;};

void Wcs::initWcsLib ( geom::Point2D const &  crval, geom::Point2D const &  crpix, Eigen::Matrix2d const &  CD, std::string const &  ctype1, std::string const &  ctype2, double  equinox, std::string const &  raDecSys, std::string const &  cunits1, std::string const &  cunits2  )

private

Manually initialise a wcs struct using values passed by the constructor.

Parameters

crval	The sky position of the reference point
crpix	The pixel position corresponding to crval in LSST units
CD	Matrix describing transformations from pixel to sky positions
ctype1	Projection system used (see description of Wcs)
ctype2	Projection system used (see description of Wcs)
equinox	Equinox of coordinate system, eg 2000 (Julian) or 1950 (Besselian)
raDecSys	System used to describe right ascension or declination, e.g FK4, FK5 or ICRS
cunits1	Units of sky position. One of deg, arcmin or arcsec
cunits2	Units of sky position. One of deg, arcmin or arcsec

Definition at line 371 of file Wcs.cc.

                                                                           {

    //Check CD is a valid size
    if( (CD.rows() != 2) || (CD.cols() != 2) ) {
        string msg = "CD is not a 2x2 matrix";
        throw LSST_EXCEPT(except::InvalidParameterError, msg);
    }

    //Check that cunits are legitimate values
    bool isValid = (cunits1 == "deg");
    isValid |= (cunits1 == "arcmin");
    isValid |= (cunits1 == "arcsec");
    isValid |= (cunits1 == "mas");
    
    if (!isValid) {
        string msg =  "CUNITS1 must be one of {deg|arcmin|arcsec|mas}";
        throw LSST_EXCEPT(except::InvalidParameterError, msg);
    }        

    isValid = (cunits2 == "deg");
    isValid |= (cunits2 == "arcmin");
    isValid |= (cunits2 == "arcsec");
    isValid |= (cunits2 == "mas");
    
    if (!isValid) {
        string msg =  "CUNITS2 must be one of {deg|arcmin|arcsec|mas}";
        throw LSST_EXCEPT(except::InvalidParameterError, msg);
    }        

    //Initialise the wcs struct
    _wcsInfo = static_cast<struct wcsprm *>(malloc(sizeof(struct wcsprm)));
    if (_wcsInfo == NULL) {
        throw LSST_EXCEPT(except::MemoryError, "Cannot allocate WCS info");
    }

    _wcsInfo->flag = -1;
    int status = wcsini(true, 2, _wcsInfo);   //2 indicates a naxis==2, a two dimensional image
    if(status != 0) {
        throw LSST_EXCEPT(except::MemoryError,
                          (boost::format("Failed to allocate memory with wcsini. Status %d: %s") %
                           status % wcs_errmsg[status] ).str());
    }
    
    //Set crval, crpix and CD. Internally to the class, we use fits units for consistency with
    //wcslib.
    _wcsInfo->crval[0] = crval.getX();
    _wcsInfo->crval[1] = crval.getY();
    _wcsInfo->crpix[0] = crpix.getX() + lsstToFitsPixels;
    _wcsInfo->crpix[1] = crpix.getY() + lsstToFitsPixels;

    //Set the CD matrix
    for (int i=0; i<2; ++i) {
        for (int j=0; j<2; ++j) {
            _wcsInfo->cd[(2*i) + j] = CD(i,j);
        }
    }

    //Specify that we have a CD matrix, but no PC or CROTA
    _wcsInfo->altlin = 2;
    _wcsInfo->flag   = 0;   //values have been updated

    //This is a work around for what I think is a bug in wcslib. ->types is neither
    //initialised or set to NULL by default, so if I try to delete a Wcs object,
    //wcslib then attempts to free non-existent space, and the code can crash.
    _wcsInfo->types = NULL;

    //Set the coordinate system
    strncpy(_wcsInfo->ctype[0], ctype1.c_str(), STRLEN);
    strncpy(_wcsInfo->ctype[1], ctype2.c_str(), STRLEN);
    strncpy(_wcsInfo->radesys, raDecSys.c_str(), STRLEN);
    _wcsInfo->equinox = equinox;
    
    //Set the units
    strncpy(_wcsInfo->cunit[0], cunits1.c_str(), STRLEN);
    strncpy(_wcsInfo->cunit[1], cunits2.c_str(), STRLEN);
    
    _nWcsInfo = 1;   //Specify that we have only one coordinate representation

    //Tell wcslib that we are need to set up internal values
    status=wcsset(_wcsInfo);
    if(status != 0) {
        throw LSST_EXCEPT(except::RuntimeError,
                          (boost::format("Failed to setup wcs structure with wcsset. Status %d: %s") %
                           status % wcs_errmsg[status] ).str());

    }
}

void Wcs::initWcsLibFromFits ( boost::shared_ptr< lsst::daf::base::PropertySet const > const &  fitsMetadata )

protected

Parse a fits header, extract the relevant metadata and create a Wcs object.

Access control for the input header

We want to hack up the input, and in order to do so we need to do a deep copy on it. We only want to do that copy once, and would like to avoid doing it altogether.

Return a readable version of the metadata

Return a writable version of the metadata

Ctor

Definition at line 172 of file Wcs.cc.

                                                                           {
    class HeaderAccess {
    public:
        CONST_PTR(lsst::daf::base::PropertySet) const& toRead() { return _constHeader; }
        PTR(lsst::daf::base::PropertySet) const& toWrite() {
            if (!_hackHeader) {
                _hackHeader = _constHeader->deepCopy();
                _constHeader = _hackHeader;
            }
            return _hackHeader;
        }
        
        HeaderAccess(CONST_PTR(lsst::daf::base::PropertySet) const& header) :
            _constHeader(header), _hackHeader() {}

    private:
        CONST_PTR(lsst::daf::base::PropertySet) _constHeader;
        PTR(lsst::daf::base::PropertySet) _hackHeader;
    };

    HeaderAccess access(header);

    // Some headers (e.g. SDSS ones from FNAL) have EQUINOX as a string.  Fix this,
    // as wcslib 4.4.4 refuses to handle it
    {
        std::string const& key = "EQUINOX";
        if (access.toRead()->exists(key) && access.toRead()->typeOf(key) == typeid(std::string)) {
            double equinox = ::atof(access.toRead()->getAsString(key).c_str());
            access.toWrite()->set(key, equinox);
        }
    }

    //Check header isn't empty
    int nCards = lsst::afw::formatters::countFitsHeaderCards(access.toRead());
    if (nCards <= 0) {
        string msg = "Could not parse FITS WCS: no header cards found";
        throw LSST_EXCEPT(except::InvalidParameterError, msg);
    }

    // Scamp produces PVi_xx header cards that are inconsistent with WCS Paper 2
    // and cause WCSLib to choke.  Aggressively, rename all PV keywords to X_PV
    for (int j=1; j<3; j++) {
        for (int i=0; i<=99; i++) {
            std::string key = (boost::format("PV%i_%i") % j % i).str();
            if (!access.toRead()->exists(key)) {
                break;
            }
            double val = access.toRead()->getAsDouble(key);
            access.toWrite()->remove(key);
            access.toWrite()->add("X_"+key, val);
        }
    }

    //While the standard does not insist on CRVAL and CRPIX being present, it 
    //is almost certain their absence indicates a problem.   
    //Check for CRPIX
    if( !access.toRead()->exists("CRPIX1") && !access.toRead()->exists("CRPIX1a")) {
        string msg = "Neither CRPIX1 not CRPIX1a found";
        throw LSST_EXCEPT(except::InvalidParameterError, msg);
    }

    if( !access.toRead()->exists("CRPIX2") && !access.toRead()->exists("CRPIX2a")) {
        string msg = "Neither CRPIX2 not CRPIX2a found";
        throw LSST_EXCEPT(except::InvalidParameterError, msg);
    }

    //And the same for CRVAL
    if( !access.toRead()->exists("CRVAL1") && !access.toRead()->exists("CRVAL1a")) {
        string msg = "Neither CRVAL1 not CRVAL1a found";
        throw LSST_EXCEPT(except::InvalidParameterError, msg);
    }

    if( !access.toRead()->exists("CRVAL2") && !access.toRead()->exists("CRVAL2a")) {
        string msg = "Neither CRVAL2 not CRVAL2a found";
        throw LSST_EXCEPT(except::InvalidParameterError, msg);
    }
    /*
     * According to Greisen and Calabretta (A&A 395, 1061–1075 (2002)) it's illegal to mix PCi_j and CDi_j
     * headers; unfortunately Subaru puts both in its headers.  It actually uses PC001002 instead of PC1_2
     * (dating to a proposed FITS standard from 1996) and at least sometimes fails to include CDELT[12],
     * so the CD and PC matrices are inconsistent
     *
     * If we detect any part of a CD matrix, delete all PC matrices
     */
    if(access.toRead()->exists("CD1_1") || access.toRead()->exists("CD1_2") ||
       access.toRead()->exists("CD2_1") || access.toRead()->exists("CD2_2")) {
        for (int i = 1; i <= 2; ++i) {
            for (int j = 1; j <= 2; ++j) {
                std::string key = (boost::format("PC%i_%i") % j % i).str();
                if (access.toRead()->exists(key)) {
                    double const val = access.toRead()->getAsDouble(key);
                    access.toWrite()->remove(key);
                    access.toWrite()->add("X_" + key, val);
                }

                key = (boost::format("PC%03d%03d") % j % i).str();
                if (access.toRead()->exists(key)) {
                    double const val = access.toRead()->getAsDouble(key);
                    access.toWrite()->remove(key);
                    access.toWrite()->add("X_" + key, val);
                }
            }
        }
    }

    //Pass the header into wcslib's formatter to extract & setup the Wcs. First need
    //to convert to a C style string, so the compile doesn't complain about constness
    std::string metadataStr = lsst::afw::formatters::formatFitsProperties(access.toRead());
    // We own the data, and wcslib is slack about constness, so no qualms with casting away const
    char *hdrString = const_cast<char*>(metadataStr.c_str());
    //printf("wcspih string:\n%s\n", hdrString);
    
    nCards = lsst::afw::formatters::countFitsHeaderCards(access.toRead()); // we may have dropped some
    int pihStatus = wcspih(hdrString, nCards, _relax, _wcshdrCtrl, &_nReject, &_nWcsInfo, &_wcsInfo);

    if (pihStatus != 0) {
        throw LSST_EXCEPT(except::RuntimeError,
                          (boost::format("Could not parse FITS WCS: wcspih status = %d (%s)") %
                           pihStatus % wcs_errmsg[pihStatus]).str());
    }    

    //Run wcsfix on _wcsInfo to try and fix any problems it knows about.
    const int *naxes = NULL;            // should be {NAXIS1, NAXIS2, ...} to check cylindrical projections
    int stats[NWCSFIX];                 // status returns from wcsfix
    int fixStatus = wcsfix(_wcsfixCtrl, naxes, _wcsInfo, stats);
    if (fixStatus != 0) {
        std::stringstream errStream;
        errStream << "Could not parse FITS WCS: wcsfix failed " << std::endl;
        for (int ii = 0; ii < NWCSFIX; ++ii) {
            if (stats[ii] >= 0) {
                errStream << "\t" << ii << ": " << stats[ii] << " " << wcsfix_errmsg[stats[ii]] << std::endl;
            } else {
                errStream << "\t" << ii << ": " << stats[ii] << std::endl;
            }
        }
    }
    
    //The Wcs standard requires a default value for RADESYS if the keyword
    //doesn't exist in header, but wcslib doesn't set it. So we do so here. This code 
    //conforms to Calabretta & Greisen 2002 \S 3.1
    if (!(access.toRead()->exists("RADESYS") || access.toRead()->exists("RADESYSa"))) {
        // If RADECSYS exists, use that (counter to Calabretta & Greisen 2002 \S 3.1, but commonly used).
        // If equinox exist and < 1984, use FK4. If >= 1984, use FK5
        if (access.toRead()->exists("RADECSYS")) {
            strncpy(_wcsInfo->radesys, access.toRead()->getAsString("RADECSYS").c_str(), STRLEN);
        } else if (access.toRead()->exists("EQUINOX") || access.toRead()->exists("EQUINOXa")) {
            std::string const EQUINOX = access.toRead()->exists("EQUINOX") ? "EQUINOX" : "EQUINOXa";
            double const equinox = access.toRead()->getAsDouble(EQUINOX);
            if(equinox < 1984) {
                strncpy(_wcsInfo->radesys, "FK4", STRLEN);
            } else {
                strncpy(_wcsInfo->radesys, "FK5", STRLEN);
            }
        } else {
            //If Equinox doesn't exist, default to ICRS
            strncpy(_wcsInfo->radesys, "ICRS", STRLEN);
        }
    }
    // strip trailing whitespace
    {
        for(int i = strlen(_wcsInfo->radesys) - 1; i >= 0; i--) {
            if (isspace(_wcsInfo->radesys[i])) {
                _wcsInfo->radesys[i] = '\0';
            }
        }
    }
    //
    // If there are no CDi_j cards in the header, set CDi_j from PCi_j
    // CDi_j == CDELTi*PCi_j
    //
    if ((_wcsInfo->altlin & 2) == 0) {  // no CDi_j cards were found in the header
        double const *cdelt = _wcsInfo->cdelt;
        double const *pc = _wcsInfo->pc;
        double *cd = _wcsInfo->cd;

        cd[0] = cdelt[0]*pc[0];         // 1_1
        cd[1] = cdelt[0]*pc[1];         // 1_2
        cd[2] = cdelt[1]*pc[2];         // 2_1
        cd[3] = cdelt[1]*pc[3];         // 2_2
    }
}

bool Wcs::isFlipped

(

)

const

Does the Wcs follow the convention of North=Up, East=Left?

This actually just measures the sign of the determinant of the CD matrix to determine the "handedness" of the coordinate system.

Returns the orientation of the Wcs

The conventional sense for a Wcs image is to have North up and East to the left, or at least to be able to rotate the image to that orientation. It is possible to create a "flipped" Wcs, where East points right when the image is rotated such that North is up. Flipping a Wcs is akin to producing a mirror image. This function tests whether the image is flipped or not.

It does so by calculating the determinant of the CD (i.e the rotation and scaling) matrix. If this determinant is positive, then the image can be rotated to a position where increasing the right ascension and declination increases the horizontal and vertical pixel position. In this case the image is flipped.

Definition at line 696 of file Wcs.cc.

                          {
    assert(_wcsInfo);
    double det = (_wcsInfo->cd[0] * _wcsInfo->cd[3]) - (_wcsInfo->cd[1] * _wcsInfo->cd[2]);

    if (det == 0) {
        throw(LSST_EXCEPT(except::RuntimeError, "Wcs CD matrix is singular"));
    }

    return (det > 0);
}

bool lsst::afw::image::Wcs::isPersistable ( ) const

virtual

Whether the Wcs is persistable using afw::table::io archives.

Reimplemented from lsst::afw::table::io::Persistable.

Definition at line 1146 of file Wcs.cc.

                              {
    if (_wcsInfo[0].naxis != 2) return false;
    if (std::strcmp(_wcsInfo[0].cunit[0], "deg") != 0) return false;
    if (std::strcmp(_wcsInfo[0].cunit[1], "deg") != 0) return false;
    return true;
}

lsst::afw::geom::AffineTransform Wcs::linearizePixelToSky	(	coord::Coord const &	coord,
		geom::AngleUnit	skyUnit = geom::degrees
	)		const

Return the local linear approximation to Wcs::pixelToSky at a point given in sky coordinates.

The local linear approximation is defined such the following is true (ignoring floating-point errors):

* wcs.linearizePixelToSky(sky, skyUnit)(wcs.skyToPixel(sky)) == sky.getPosition(skyUnit);
* 

(recall that AffineTransform::operator() is matrix multiplication with the augmented point (x,y,1)).

This is currently implemented as a numerical derivative, but we should specialise the Wcs class (or rather its implementation) to handle "simple" cases such as TAN-SIP analytically

Parameters

[in]	coord	Position in sky coordinates where transform is desired.
[in]	skyUnit	Units to use for sky coordinates; units of matrix elements will be skyUnits/pixel.

Definition at line 977 of file Wcs.cc.

        {
    return linearizePixelToSkyInternal(skyToPixel(coord), coord, skyUnit);
}

lsst::afw::geom::AffineTransform Wcs::linearizePixelToSky	(	geom::Point2D const &	pix,
		geom::AngleUnit	skyUnit = geom::degrees
	)		const

Return the local linear approximation to Wcs::pixelToSky at a point given in pixel coordinates.

The local linear approximation is defined such the following is true (ignoring floating-point errors):

* wcs.linearizePixelToSky(pix, skyUnit)(pix) == wcs.pixelToSky(pix).getPosition(skyUnit)
* 

(recall that AffineTransform::operator() is matrix multiplication with the augmented point (x,y,1)).

This is currently implemented as a numerical derivative, but we should specialise the Wcs class (or rather its implementation) to handle "simple" cases such as TAN-SIP analytically

Parameters

[in]	pix	Position in pixel coordinates where transform is desired.
[in]	skyUnit	Units to use for sky coordinates; units of matrix elements will be skyUnits/pixel.

Definition at line 983 of file Wcs.cc.

        {
    return linearizePixelToSkyInternal(pix, *pixelToSky(pix), skyUnit);
}

lsst::afw::geom::AffineTransform Wcs::linearizePixelToSkyInternal ( geom::Point2D const &  pix, coord::Coord const &  coord, geom::AngleUnit  skyUnit  ) const

protectedvirtual

Definition at line 994 of file Wcs.cc.

        {
    //
    // Figure out the (0, 0), (0, 1), and (1, 0) ra/dec coordinates of the corners of a square drawn in pixel
    // It'd be better to centre the square at sky00, but that would involve another conversion between sky and
    // pixel coordinates so I didn't bother
    //
    const double side = 10;             // length of the square's sides in pixels
    GeomPoint const sky00 = coord.getPosition(skyUnit);
    typedef std::pair<lsst::afw::geom::Angle, lsst::afw::geom::Angle> AngleAngle;
    AngleAngle const dsky10 = coord.getTangentPlaneOffset(*pixelToSky(pix00 + afwGeom::Extent2D(side, 0)));
    AngleAngle const dsky01 = coord.getTangentPlaneOffset(*pixelToSky(pix00 + afwGeom::Extent2D(0, side)));

    Eigen::Matrix2d m;
    m(0, 0) = dsky10.first.asAngularUnits(skyUnit)/side;
    m(0, 1) = dsky01.first.asAngularUnits(skyUnit)/side;
    m(1, 0) = dsky10.second.asAngularUnits(skyUnit)/side;
    m(1, 1) = dsky01.second.asAngularUnits(skyUnit)/side;

    Eigen::Vector2d sky00v;
    sky00v << sky00.getX(), sky00.getY();
    Eigen::Vector2d pix00v;
    pix00v << pix00.getX(), pix00.getY();
    //return lsst::afw::geom::AffineTransform(m, lsst::afw::geom::Extent2D(sky00v - m * pix00v));
    return lsst::afw::geom::AffineTransform(m, (sky00v - m * pix00v));
}

lsst::afw::geom::AffineTransform Wcs::linearizeSkyToPixel	(	coord::Coord const &	coord,
		geom::AngleUnit	skyUnit = geom::degrees
	)		const

Return the local linear approximation to Wcs::skyToPixel at a point given in sky coordinates.

The local linear approximation is defined such the following is true (ignoring floating-point errors):

* wcs.linearizeSkyToPixel(sky, skyUnit)(sky.getPosition(skyUnit)) == wcs.skyToPixel(sky)
* 

(recall that AffineTransform::operator() is matrix multiplication with the augmented point (x,y,1)).

This is currently implemented as a numerical derivative, but we should specialise the Wcs class (or rather its implementation) to handle "simple" cases such as TAN-SIP analytically

Parameters

[in]	coord	Position in sky coordinates where transform is desired.
[in]	skyUnit	Units to use for sky coordinates; units of matrix elements will be pixels/skyUnit.

Definition at line 1024 of file Wcs.cc.

        {
    return linearizeSkyToPixelInternal(skyToPixel(coord), coord, skyUnit);
}

lsst::afw::geom::AffineTransform Wcs::linearizeSkyToPixel	(	geom::Point2D const &	pix,
		geom::AngleUnit	skyUnit = geom::degrees
	)		const

Return the local linear approximation to Wcs::skyToPixel at a point given in pixel coordinates.

The local linear approximation is defined such the following is true (ignoring floating-point errors):

* wcs.linearizeSkyToPixel(pix, skyUnit)(wcs.pixelToSky(pix).getPosition(skyUnit)) == pix
* 

(recall that AffineTransform::operator() is matrix multiplication with the augmented point (x,y,1)).

This is currently implemented as a numerical derivative, but we should specialise the Wcs class (or rather its implementation) to handle "simple" cases such as TAN-SIP analytically

Parameters

[in]	pix	Position in pixel coordinates where transform is desired.
[in]	skyUnit	Units to use for sky coordinates; units of matrix elements will be pixels/skyUnit.

Definition at line 1031 of file Wcs.cc.

        {
    return linearizeSkyToPixelInternal(pix, *pixelToSky(pix), skyUnit);
}

lsst::afw::geom::AffineTransform Wcs::linearizeSkyToPixelInternal ( geom::Point2D const &  pix, coord::Coord const &  coord, geom::AngleUnit  skyUnit  ) const

protectedvirtual

Implementation for the overloaded public linearizeSkyToPixel methods, requiring both a pixel coordinate and the corresponding sky coordinate.

Definition at line 1042 of file Wcs.cc.

        {
    lsst::afw::geom::AffineTransform inverse = linearizePixelToSkyInternal(pix00, coord, skyUnit);
    return inverse.invert();
}

CoordPtr Wcs::makeCorrectCoord ( geom::Angle  sky0, geom::Angle  sky1  ) const

protected

Given a sky position, use the values stored in ctype and radesys to return the correct sub-class of Coord.

Definition at line 922 of file Wcs.cc.

                                                                                 {

    //Construct a coord object of the correct type
    int const ncompare = 4;                       // we only care about type's first 4 chars
    char *type = _wcsInfo->ctype[0];
    char *radesys = _wcsInfo->radesys;
    double equinox = _wcsInfo->equinox;

    if (strncmp(type, "RA--", ncompare) == 0) { // Our default.  If it's often something else, consider
        ;                                       // using an tr1::unordered_map
        if(strcmp(radesys, "ICRS") == 0) {
            return afwCoord::makeCoord(afwCoord::ICRS, sky0, sky1);
        }
        if(strcmp(radesys, "FK5") == 0) {
            return afwCoord::makeCoord(afwCoord::FK5, sky0, sky1, equinox);
        } else {   
            throw LSST_EXCEPT(except::RuntimeError,
                              (boost::format("Can't create Coord object: Unrecognised radesys %s") %
                               radesys).str());
        }

    } else if (strncmp(type, "GLON", ncompare) == 0) {
        return afwCoord::makeCoord(afwCoord::GALACTIC, sky0, sky1);   
    } else if (strncmp(type, "ELON", ncompare) == 0) {
        return afwCoord::makeCoord(afwCoord::ECLIPTIC, sky0, sky1, equinox);
    } else if (strncmp(type, "DEC-", ncompare) == 0) {
        //check for the case where the ctypes are swapped. Note how sky0 and sky1 are swapped as well

        //Our default
        if(strcmp(radesys, "ICRS") == 0) {
            return afwCoord::makeCoord(afwCoord::ICRS, sky1, sky0);
        }
        if(strcmp(radesys, "FK5") == 0) {
            return afwCoord::makeCoord(afwCoord::FK5, sky1, sky0, equinox);
        } else {   
            throw LSST_EXCEPT(except::RuntimeError,
                              (boost::format("Can't create Coord object: Unrecognised radesys %s") %
                               radesys).str());
        }
    } else if (strncmp(type, "GLAT", ncompare) == 0) {
        return afwCoord::makeCoord(afwCoord::GALACTIC, sky1, sky0);   
    } else if (strncmp(type, "ELAT", ncompare) == 0) {
        return afwCoord::makeCoord(afwCoord::ECLIPTIC, sky1, sky0, equinox);
    } else {
    //Give up in disgust
        throw LSST_EXCEPT(except::RuntimeError,
                          (boost::format("Can't create Coord object: Unrecognised sys %s") %
                           type).str());
    }
    
    //Can't get here
    assert(0);
}

bool lsst::afw::image::Wcs::operator!= ( Wcs const &  other ) const

inline

Definition at line 135 of file Wcs.h.

{ return !(*this == other); }

Wcs& lsst::afw::image::Wcs::operator= ( const Wcs &  )

protected

bool Wcs::operator==

(

Wcs const &

other

)

const

Definition at line 496 of file Wcs.cc.

                                            {
    if (&other == this) return true;
    // We do a bidirectional test with a virtual member function in case one of us is a derived
    // class with members we don't know about here.
    // This is not the most efficient possible implementation, but I think it's the easiest one
    // with which to ensure correctness, and I think that's more important in this case.
    return this->_isSubset(other) && other._isSubset(*this);
}

double Wcs::pixArea

(

lsst::afw::geom::Point2D

pix00

)

const

Sky area covered by a pixel at position pix00 in units of square degrees.

Sky area covered by a pixel at position pix00. In units of square degrees.

Parameters

pix00

The pixel point where the area is desired

Definition at line 712 of file Wcs.cc.

                           {
    //
    // Figure out the (0, 0), (0, 1), and (1, 0) ra/dec coordinates of the corners of a square drawn in pixel
    // It'd be better to centre the square at sky00, but that would involve another conversion between sky and
    // pixel coordinates so I didn't bother
    //
    const double side = 1;             // length of the square's sides in pixels

    // Work in 3-space to avoid RA wrapping and pole issues.
    afwGeom::Point3D v0 = pixelToSky(pix00)->getVector();

    // Step by "side" in x and y pixel directions...
    GeomPoint px(pix00);
    GeomPoint py(pix00);
    px.shift(afwGeom::Extent2D(side, 0));
    py.shift(afwGeom::Extent2D(0, side));
    // Push the points through the WCS, and find difference in 3-space.
    afwGeom::Extent3D dx = pixelToSky(px)->getVector() - v0;
    afwGeom::Extent3D dy = pixelToSky(py)->getVector() - v0;

    // Compute |cross product| = area of parallelogram with sides dx,dy
    // FIXME -- this is slightly incorrect -- it's making the small-angle
    // approximation, taking the distance *through* the unit sphere
    // rather than over its surface.
    // This is in units of ~radians^2
    double area = sqrt(square(dx[1]*dy[2] - dx[2]*dy[1]) +
                       square(dx[2]*dy[0] - dx[0]*dy[2]) +
                       square(dx[0]*dy[1] - dx[1]*dy[0]));

    return area / square(side) * square(180. / afwGeom::PI);
}

afwGeom::Angle Wcs::pixelScale

(

)

const

Returns the pixel scale [Angle/pixel].

Definition at line 745 of file Wcs.cc.

                                   {
    return sqrt(pixArea(getPixelOrigin())) * afwGeom::degrees;
}

CoordPtr Wcs::pixelToSky	(	double	pixel1,
		double	pixel2
	)		const

Convert from celestial coordinates to pixel coordinates.

Convert from pixel position to sky coordinates (e.g ra/dec)

Convert a pixel position (e.g x,y) to a celestial coordinate (e.g ra/dec). The output coordinate system depends on the values of CTYPE used to construct the object. For ra/dec, the CTYPES should be RA–TAN and DEC-TAN.

Definition at line 894 of file Wcs.cc.

                                                           {
    assert(_wcsInfo);

    afwGeom::Angle skyTmp[2];
    pixelToSkyImpl(pixel1, pixel2, skyTmp);
    return makeCorrectCoord(skyTmp[0], skyTmp[1]);
}

CoordPtr Wcs::pixelToSky

(

lsst::afw::geom::Point2D const &

pixel

)

const

Convert from celestial coordiantes to pixel coordinates.

Convert from pixel position to sky coordinates (e.g ra/dec)

Convert a pixel position (e.g x,y) to a celestial coordinate (e.g ra/dec). The output coordinate system depends on the values of CTYPE used to construct the object. For ra/dec, the CTYPES should be RA–TAN and DEC-TAN.

Definition at line 885 of file Wcs.cc.

                                                      {
    return pixelToSky(pixel.getX(), pixel.getY());
}

void Wcs::pixelToSky	(	double	pixel1,
		double	pixel2,
		geom::Angle &	sky1,
		geom::Angle &	sky2
	)		const

Convert from celestial coordiantes to pixel coordinates.

Convert from pixel position to sky coordinates (e.g ra/dec)

Note

This routine is designed for the knowledgeable user in need of performance; it's safer to call the version that returns a PTR(Coord).

Convert a pixel position (e.g x,y) to a celestial coordinate (e.g ra/dec)

Note

This routine is designed for the knowledgeable user in need of performance; it's safer to call the version that returns a CoordPtr

Definition at line 909 of file Wcs.cc.

                                                                                               {
    afwGeom::Angle skyTmp[2];
    // HACK -- we shouldn't need to initialize these -- pixelToSkyImpl() sets them unless an
    // exception is thrown -- but be safe.
    skyTmp[0] = 0. * afwGeom::radians;
    skyTmp[1] = 0. * afwGeom::radians;
    pixelToSkyImpl(pixel1, pixel2, skyTmp);
    sky1 = skyTmp[0];
    sky2 = skyTmp[1];
}

void Wcs::pixelToSkyImpl ( double  pixel1, double  pixel2, geom::Angle  skyTmp[2]  ) const

privatevirtual

Reimplemented in lsst::afw::image::TanWcs.

Definition at line 857 of file Wcs.cc.

{
    assert(_wcsInfo);
    
    // wcslib assumes 1-indexed coordinates
    double pixTmp[2] = { pixel1 - lsst::afw::image::PixelZeroPos + lsstToFitsPixels,
                         pixel2 - lsst::afw::image::PixelZeroPos + lsstToFitsPixels}; 
    double imgcrd[2];
    double phi, theta;
    
    double sky[2];
    int status = 0;
    status = wcsp2s(_wcsInfo, 1, 2, pixTmp, imgcrd, &phi, &theta, sky, &status);
    if (status > 0) {
        throw LSST_EXCEPT(except::RuntimeError,
            (boost::format("Error: wcslib returned a status code of %d at pixel %s, %s: %s") %
            status % pixel1 % pixel2 % wcs_errmsg[status]).str());
    }
    // FIXME -- _wcsInfo.lat, _wcsInfo.lng ?
    skyTmp[0] = sky[0] * afwGeom::degrees;
    skyTmp[1] = sky[1] * afwGeom::degrees;
}

void Wcs::rotateImageBy90 ( int  nQuarter, lsst::afw::geom::Extent2I  dimensions  ) const

virtual

Definition at line 627 of file Wcs.cc.

                                                                        {
    assert(_wcsInfo);

    while (nQuarter < 0 ) {
        nQuarter += 4;
    }

    
    int const naxis = _wcsInfo->naxis;

    //If naxis != 2, I'm not sure if any of what follows is correct
    assert(naxis == 2);
    double a = _wcsInfo->cd[0];
    double b = _wcsInfo->cd[1];
    double c = _wcsInfo->cd[2];
    double d = _wcsInfo->cd[3];
    double crpx = _wcsInfo->crpix[0];
    double crpy = _wcsInfo->crpix[1];
    switch (nQuarter%4) {
        case 0:
            break;
        case 1:
            _wcsInfo->cd[0] = -b;
            _wcsInfo->cd[1] = a;
            _wcsInfo->cd[2] = -d;
            _wcsInfo->cd[3] = c;
            _wcsInfo->crpix[0] = -crpy + dimensions.getY();
            _wcsInfo->crpix[1] = crpx;
            break;
        case 2:
            _wcsInfo->cd[0] = -a;
            _wcsInfo->cd[1] = -b;
            _wcsInfo->cd[2] = -c;
            _wcsInfo->cd[3] = -d;
            _wcsInfo->crpix[0] = -crpx + dimensions.getX();
            _wcsInfo->crpix[1] = -crpy + dimensions.getY();
            break;
        case 3:
            _wcsInfo->cd[0] = b;
            _wcsInfo->cd[1] = -a;
            _wcsInfo->cd[2] = d;
            _wcsInfo->cd[3] = -c;
            _wcsInfo->crpix[0] = crpy;
            _wcsInfo->crpix[1] = -crpx + dimensions.getX();
            break;
    }

    // tells libwcs to invalidate cached data, since transformation has been modified
    _wcsInfo->flag = 0;
}

void lsst::afw::image::Wcs::shiftReferencePixel ( geom::Extent2D const &  d )

inline

Definition at line 283 of file Wcs.h.

{shiftReferencePixel(d.getX(), d.getY());}

void Wcs::shiftReferencePixel	(	double	dx,
		double	dy
	)

Move the pixel reference position by (dx, dy) Used when persisting and retrieving sub-images. The lsst convention is that Wcs returns pixel position (which is based on position in the parent image), but the fits convention is to return pixel index (which is bases on position in the sub-image). In order that the fits files we create make sense to other fits viewers, we change to the fits convention when writing out images.

Definition at line 1205 of file Wcs.cc.

                                                  {
    assert(_wcsInfo);
    _wcsInfo->crpix[0] += dx;
    _wcsInfo->crpix[1] += dy;

    // tells libwcs to invalidate cached data, since transformation has been modified
    _wcsInfo->flag = 0;
}

GeomPoint Wcs::skyToIntermediateWorldCoord

(

coord::Coord const &

coord

)

const

Convert from sky coordinates (e.g ra/dec) to intermediate world coordinates.

Intermediate world coordinates are in DEGREES.

Definition at line 812 of file Wcs.cc.

                                                                                {
    assert(_wcsInfo);

    afwCoord::Coord::Ptr sky = convertCoordToSky(coord);
    double skyTmp[2];
    double imgcrd[2];
    double phi, theta;
    double pixTmp[2];

    /*
     printf("skyToIWC: _coordSystem = %i\n", (int)_coordSystem);
     printf("coord (%.3f, %.3f)\n", coord->getLongitude().asDegrees(), coord->getLatitude().asDegrees());
     printf("->sky (%.3f, %.3f)\n", sky->getLongitude().asDegrees(), sky->getLatitude().asDegrees());
     */

    skyTmp[_wcsInfo->lng] = sky->getLongitude().asDegrees();
    skyTmp[_wcsInfo->lat] = sky->getLatitude() .asDegrees();

    //Estimate pixel coordinates
    int stat[1];
    int status = 0;
    imgcrd[0] = imgcrd[1] = -1e6;
    /*
     printf("  skyTmp[] = (%.3f, %.3f)\n", skyTmp[0], skyTmp[1]);
     printf("  _wcsInfo->lng,lat = %i, %i\n", _wcsInfo->lng, _wcsInfo->lat);
     */
    status = wcss2p(_wcsInfo, 1, 2, skyTmp, &phi, &theta, imgcrd, pixTmp, stat);
    if (status > 0) {
        throw LSST_EXCEPT(except::RuntimeError,
            (boost::format("Error: wcslib returned a status code of %d at sky %s, %s deg: %s") %
            status % skyTmp[0] % skyTmp[1] % wcs_errmsg[status]).str());
    }
    /*
     printf("->iwc (%.3f, %.3f)\n", imgcrd[0], imgcrd[1]);
     printf("-> pix (%.2f, %.2f)\n", pixTmp[0], pixTmp[1]);
     afwCoord::Coord::Ptr crval = getSkyOrigin();
     printf("(crval is (%.3f, %.3f))\n", crval->getLongitude().asDegrees(), crval->getLatitude().asDegrees());
     */
    return GeomPoint(imgcrd[0], imgcrd[1]); 
}

GeomPoint Wcs::skyToPixel	(	geom::Angle	sky1,
		geom::Angle	sky2
	)		const

Convert from sky coordinates (e.g ra/dec) to pixel positions.

ASSUMES the angles are in the appropriate coordinate system for this WCS.

Convert a sky position (e.g RA/Dec) to a pixel position. The exact meaning of sky1, sky2 and the return value depend on the properties of the wcs (i.e the values of CTYPE1 and CTYPE2), but the inputs are usually RA/Dec. The outputs are x and y pixel position.

Definition at line 808 of file Wcs.cc.

                                                                    {
    return skyToPixelImpl(sky1, sky2);
}

GeomPoint Wcs::skyToPixel

(

coord::Coord const &

coord

)

const

Convert from sky coordinates (e.g ra/dec) to pixel positions.

Definition at line 791 of file Wcs.cc.

                                                               {
    afwCoord::Coord::Ptr sky = convertCoordToSky(coord);
    return skyToPixelImpl(sky->getLongitude(), sky->getLatitude());
}

GeomPoint Wcs::skyToPixelImpl ( geom::Angle  sky1, geom::Angle  sky2  ) const

privatevirtual

Reimplemented in lsst::afw::image::TanWcs.

Definition at line 752 of file Wcs.cc.

                                     {
    assert(_wcsInfo);

    double skyTmp[2];
    double imgcrd[2];
    double phi, theta;
    double pixTmp[2];
    /*
     printf("_skyCoordsReversed: %c\n", (_skyCoordsReversed ? 'T' : 'F'));
     printf("wcsinfo.lat: %i,  lng: %i\n", _wcsInfo->lat, _wcsInfo->lng);
     */
    // WCSLib is smart enough to notice and handle crazy SDSS CTYPE1 = DEC--TAN,
    // by recording the indices of the long and lat coordinates.
    skyTmp[_wcsInfo->lng] = sky1.asDegrees();
    skyTmp[_wcsInfo->lat] = sky2.asDegrees();

    int stat[1];
    int status = 0;
    status = wcss2p(_wcsInfo, 1, 2, skyTmp, &phi, &theta, imgcrd, pixTmp, stat);
    if (status == 9) {
        throw LSST_EXCEPT(except::DomainError,
            (boost::format("sky coordinates %s, %s degrees is not valid for this WCS")
             % sky1.asDegrees() % sky2.asDegrees()
             ).str()
        );
    }
    if (status > 0) {
        throw LSST_EXCEPT(except::RuntimeError,
            (boost::format("Error: wcslib returned a status code of %d at sky %s, %s deg: %s") %
            status % sky1.asDegrees() % sky2.asDegrees() % wcs_errmsg[status]).str());
    }

    // wcslib assumes 1-indexed coords
    return afwGeom::Point2D(pixTmp[0] + lsst::afw::image::PixelZeroPos + fitsToLsstPixels,
                            pixTmp[1] + lsst::afw::image::PixelZeroPos + fitsToLsstPixels); 
}

void lsst::afw::image::Wcs::write ( OutputArchiveHandle &  handle ) const

protectedvirtual

Write the object to one or more catalogs.

The handle object passed to this function provides an interface for adding new catalogs and adding nested objects to the same archive (while checking for duplicates). See OutputArchiveHandle for more information.

Reimplemented from lsst::afw::table::io::Persistable.

Reimplemented in lsst::afw::image::TanWcs.

Definition at line 1128 of file Wcs.cc.

                                                  {
    WcsPersistenceHelper const & keys = WcsPersistenceHelper::get();
    afw::table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
    PTR(afw::table::BaseRecord) record = catalog.addNew();
    record->set(keys.crval, getSkyOrigin()->getPosition(afw::geom::degrees));
    record->set(keys.crpix, getPixelOrigin());
    Eigen::Matrix2d cdIn = getCDMatrix();
    Eigen::Map<Eigen::Matrix2d> cdOut((*record)[keys.cd].getData());
    cdOut = cdIn;
    record->set(keys.ctype1, std::string(_wcsInfo[0].ctype[0]));
    record->set(keys.ctype2, std::string(_wcsInfo[0].ctype[1]));
    record->set(keys.equinox, _wcsInfo[0].equinox);
    record->set(keys.radesys, std::string(_wcsInfo[0].radesys));
    record->set(keys.cunit1, std::string(_wcsInfo[0].cunit[0]));
    record->set(keys.cunit2, std::string(_wcsInfo[0].cunit[1]));
    handle.saveCatalog(catalog);
}

Friends And Related Function Documentation

Wcs::Ptr makeWcs ( boost::shared_ptr< lsst::daf::base::PropertySet > const &  fitsMetadata, bool  stripMetadata  )

friend

Create a Wcs of the correct class using a fits header.

Set stripMetadata=true to remove processed keywords from the PropertySet.

friend class WcsFactory

friend

Definition at line 305 of file Wcs.h.

Member Data Documentation

coord::CoordSystem lsst::afw::image::Wcs::_coordSystem

protected

Definition at line 359 of file Wcs.h.

int lsst::afw::image::Wcs::_nReject

protected

Definition at line 358 of file Wcs.h.

int lsst::afw::image::Wcs::_nWcsInfo

protected

Definition at line 354 of file Wcs.h.

int lsst::afw::image::Wcs::_relax

protected

Degree of permissiveness for wcspih (0 for strict); see wcshdr.h for details.

Definition at line 355 of file Wcs.h.

int lsst::afw::image::Wcs::_wcsfixCtrl

protected

Do potentially unsafe translations of non-standard unit strings? 0/1 = no/yes.

Definition at line 356 of file Wcs.h.

int lsst::afw::image::Wcs::_wcshdrCtrl

protected

Controls messages to stderr from wcshdr (0 for none); see wcshdr.h for details.

Definition at line 357 of file Wcs.h.

struct wcsprm* lsst::afw::image::Wcs::_wcsInfo

protected

Definition at line 353 of file Wcs.h.

---

The documentation for this class was generated from the following files:

• /lsst/home/lsstsw/stack/Linux64/afw/master-gae20e2b580+3/include/lsst/afw/image/Wcs.h
• /lsst/home/lsstsw/stack/Linux64/afw/master-gae20e2b580+3/src/image/Wcs.cc