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. This need not be the centre of the image. More...
lsst::afw::geom::Point2D	getPixelOrigin () const
	Returns CRPIX (corrected to LSST convention). More...
Eigen::Matrix2d	getCDMatrix () const
	Returns the CD matrix. 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
	Rotate image by nQuarter times 90 degrees. More...
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 pixel position to sky coordinates (e.g. RA/dec) More...
boost::shared_ptr< coord::Coord >	pixelToSky (lsst::afw::geom::Point2D const &pixel) const
	Convert from pixel position to sky coordinates (e.g. RA/dec) More...
void	pixelToSky (double pixel1, double pixel2, geom::Angle &sky1, geom::Angle &sky2) const
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...
virtual void	shiftReferencePixel (double dx, double dy)
	Move the pixel reference position by (dx, dy) More...
virtual void	shiftReferencePixel (geom::Extent2D const &d)
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
bool	_mayBePersistable () const
	Perform basic checks on whether *this might be persistable. More...
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 &)
virtual void	pixelToSkyImpl (double pixel1, double pixel2, geom::Angle skyTmp[2]) const
virtual geom::Point2D	skyToPixelImpl (geom::Angle sky1, geom::Angle sky2) const
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...

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 149 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 545 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 88 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 99 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 1112 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 454 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 118 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 520 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]);
}

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

protected

Perform basic checks on whether *this might be persistable.

Definition at line 1096 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;
}

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

protected

Definition at line 69 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, and lsst::afw::image::DistortedTanWcs.

Definition at line 552 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 778 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.

Reimplemented in lsst::afw::image::DistortedTanWcs.

Definition at line 591 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);
    //Origin is at (1,1).  Adjust to avoid off by one.
    if (flipLR) {
        _wcsInfo->cd[0] = -_wcsInfo->cd[0];
        _wcsInfo->cd[2] = -_wcsInfo->cd[2];
        _wcsInfo->crpix[0] = -_wcsInfo->crpix[0] + dimensions.getX() + 1;
    }
    if (flipTB) {
        _wcsInfo->cd[1] = -_wcsInfo->cd[1];
        _wcsInfo->cd[3] = -_wcsInfo->cd[3];
        _wcsInfo->crpix[1] = -_wcsInfo->crpix[1]+dimensions.getY() + 1;
    }

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

Eigen::Matrix2d Wcs::getCDMatrix

(

)

const

Returns the CD matrix.

Definition at line 573 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.

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

Definition at line 668 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.

Definition at line 1008 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 1074 of file Wcs.cc.

{ return getWcsPersistenceName(); }

GeomPoint Wcs::getPixelOrigin

(

)

const

Returns CRPIX (corrected to LSST convention).

Definition at line 565 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 1076 of file Wcs.cc.

{ return "lsst.afw.image"; }

CoordPtr Wcs::getSkyOrigin

(

)

const

Returns CRVAL. This need not be the centre of the image.

Definition at line 560 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, and lsst::afw::image::DistortedTanWcs.

Definition at line 220 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 361 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 173 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.  Furthermore, some headers (e.g. Skymapper)
    // use a string but prepend 'J':  "J2000.0" -- if we don't handle this, we deduce
    // an equinox of 0.0
    {
        std::string const& key = "EQUINOX";
        if (access.toRead()->exists(key) && access.toRead()->typeOf(key) == typeid(std::string)) {
            auto equinoxStr = access.toRead()->getAsString(key).c_str();
            double equinox = ::atof(equinoxStr[0] == 'J' ? equinoxStr + 1 : equinoxStr);
            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);
    }

    //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?

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.

Definition at line 672 of file Wcs.cc.

                          {
    // We calculate 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.
    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.

Reimplemented in lsst::afw::image::TanWcs, and lsst::afw::image::DistortedTanWcs.

Definition at line 1104 of file Wcs.cc.

                              {
    if (!_mayBePersistable()) {
        return false;
    }
    // The current table persistence only works for TAN and TAN-SIP projections
    return false;
}

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 934 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 940 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 951 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 981 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 988 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 999 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 879 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 486 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.

Parameters

pix00

The pixel point where the area is desired

Definition at line 692 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 725 of file Wcs.cc.

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

CoordPtr Wcs::pixelToSky	(	double	pix1,
		double	pix2
	)		const

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 858 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 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 854 of file Wcs.cc.

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

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

Definition at line 866 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

protectedvirtual

Reimplemented in lsst::afw::image::TanWcs, and lsst::afw::image::DistortedTanWcs.

Definition at line 831 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

Rotate image by nQuarter times 90 degrees.

Reimplemented in lsst::afw::image::DistortedTanWcs.

Definition at line 614 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];
    //Origin is at (1,1).  Adjust to avoid off by one.
    //E.g. CRPIX one pixel off the UR of the grid should go to
    //(0,0) for nQuarter=2
    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() + 1;
            _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() + 1;
            _wcsInfo->crpix[1] = -crpy + dimensions.getY() + 1;
            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() + 1;
            break;
    }

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

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

virtual

Move the pixel reference position by (dx, 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.

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.

Reimplemented in lsst::afw::image::DistortedTanWcs.

Definition at line 1164 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;
}

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

inlinevirtual

Definition at line 321 of file Wcs.h.

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

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 786 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.

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.

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

Definition at line 782 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 771 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

protectedvirtual

Reimplemented in lsst::afw::image::TanWcs, and lsst::afw::image::DistortedTanWcs.

Definition at line 732 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 1078 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 339 of file Wcs.h.

Member Data Documentation

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

protected

Definition at line 397 of file Wcs.h.

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

protected

Definition at line 396 of file Wcs.h.

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

protected

Definition at line 392 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 393 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 394 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 395 of file Wcs.h.

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

protected

Definition at line 391 of file Wcs.h.

---

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

• /home/lsstsw/stack/Linux64/afw/11.0-2-g04d2804/include/lsst/afw/image/Wcs.h
• /home/lsstsw/stack/Linux64/afw/11.0-2-g04d2804/src/image/Wcs.cc