LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
Public Types | Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Friends | List of all members
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:
lsst::daf::base::Persistable lsst::daf::base::Citizen lsst.afw.table.io::PersistableFacade< Wcs > lsst.afw.table.io::Persistable lsst.afw.image::TanWcs lsst.afw.image::DistortedTanWcs

Public Types

typedef boost::shared_ptr< WcsPtr
 
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::CoordpixelToSky (double pix1, double pix2) const
 Convert from pixel position to sky coordinates (e.g. RA/dec) More...
 
boost::shared_ptr< coord::CoordpixelToSky (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 ()
 
Citizenoperator= (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...
 
Wcsoperator= (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< WcsreadFits (fits::Fits &fitsfile)
 Read an object from an already open FITS object. More...
 
static boost::shared_ptr< WcsreadFits (std::string const &fileName, int hdu=0)
 Read an object from a regular FITS file. More...
 
static boost::shared_ptr< WcsreadFits (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

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
crvalThe sky position of the reference point
crpixThe pixel position corresponding to crval in LSST units
CDMatrix describing transformations from pixel to sky positions
ctype1Projection system used (see description of Wcs)
ctype2Projection system used (see description of Wcs)
equinoxEquinox of coordinate system, eg 2000 (Julian) or 1950 (Besselian)
raDecSysSystem used to describe right ascension or declination, e.g FK4, FK5 or ICRS
cunits1Units of sky position. One of deg, arcmin or arcsec
cunits2Units 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.

153  :
154  daf::base::Citizen(typeid(this)),
155  _wcsInfo(NULL),
156  _nWcsInfo(0),
157  _relax(0),
158  _wcsfixCtrl(0),
159  _wcshdrCtrl(0),
160  _nReject(0),
161  _coordSystem(static_cast<afwCoord::CoordSystem>(-1))
162 {
164  initWcsLib(crval, crpix, CD,
165  ctype1, ctype2,
166  equinox, raDecSys,
167  cunits1, cunits2);
168  _initWcs();
169 }
int _wcsfixCtrl
Do potentially unsafe translations of non-standard unit strings? 0/1 = no/yes.
Definition: Wcs.h:394
void _setWcslibParams()
Definition: Wcs.cc:69
coord::CoordSystem _coordSystem
Definition: Wcs.h:397
table::PointKey< double > crval
Definition: Wcs.cc:1035
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.
Definition: Wcs.cc:361
table::Key< std::string > ctype2
Definition: Wcs.cc:1039
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
table::Key< double > equinox
Definition: Wcs.cc:1040
void _initWcs()
Definition: Wcs.cc:118
int _relax
Degree of permissiveness for wcspih (0 for strict); see wcshdr.h for details.
Definition: Wcs.h:393
int _wcshdrCtrl
Controls messages to stderr from wcshdr (0 for none); see wcshdr.h for details.
Definition: Wcs.h:395
table::Key< std::string > ctype1
Definition: Wcs.cc:1038
table::PointKey< double > crpix
Definition: Wcs.cc:1036
Wcs::~Wcs ( )
virtual

Definition at line 545 of file Wcs.cc.

545  {
546  if (_wcsInfo != NULL) {
547  wcsvfree(&_nWcsInfo, &_wcsInfo);
548  }
549 }
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
lsst.afw.image::Wcs::Wcs ( )
protected

Construct an invalid Wcs given no arguments.

Definition at line 88 of file Wcs.cc.

88  :
89  daf::base::Citizen(typeid(this)),
90  _wcsInfo(NULL), _nWcsInfo(0), _relax(0), _wcsfixCtrl(0), _wcshdrCtrl(0), _nReject(0),
91  _coordSystem(static_cast<afwCoord::CoordSystem>(-1)) {
93  _initWcs();
94 }
int _wcsfixCtrl
Do potentially unsafe translations of non-standard unit strings? 0/1 = no/yes.
Definition: Wcs.h:394
void _setWcslibParams()
Definition: Wcs.cc:69
coord::CoordSystem _coordSystem
Definition: Wcs.h:397
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
void _initWcs()
Definition: Wcs.cc:118
int _relax
Degree of permissiveness for wcspih (0 for strict); see wcshdr.h for details.
Definition: Wcs.h:393
int _wcshdrCtrl
Controls messages to stderr from wcshdr (0 for none); see wcshdr.h for details.
Definition: Wcs.h:395
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.

99  :
100  daf::base::Citizen(typeid(this)),
101  _wcsInfo(NULL),
102  _nWcsInfo(0),
103  _relax(0),
104  _wcsfixCtrl(0),
105  _wcshdrCtrl(0),
106  _nReject(0),
107  _coordSystem(static_cast<afwCoord::CoordSystem>(-1))
108 {
110 
111  initWcsLibFromFits(fitsMetadata);
112  _initWcs();
113 }
int _wcsfixCtrl
Do potentially unsafe translations of non-standard unit strings? 0/1 = no/yes.
Definition: Wcs.h:394
void _setWcslibParams()
Definition: Wcs.cc:69
coord::CoordSystem _coordSystem
Definition: Wcs.h:397
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.
Definition: Wcs.cc:173
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
void _initWcs()
Definition: Wcs.cc:118
int _relax
Degree of permissiveness for wcspih (0 for strict); see wcshdr.h for details.
Definition: Wcs.h:393
int _wcshdrCtrl
Controls messages to stderr from wcshdr (0 for none); see wcshdr.h for details.
Definition: Wcs.h:395
lsst.afw.image::Wcs::Wcs ( afw::table::BaseRecord const &  record)
explicitprotected

Definition at line 1112 of file Wcs.cc.

1112  :
1113  daf::base::Citizen(typeid(this)),
1114  _wcsInfo(NULL),
1115  _nWcsInfo(0),
1116  _relax(0),
1117  _wcsfixCtrl(0),
1118  _wcshdrCtrl(0),
1119  _nReject(0),
1120  _coordSystem(static_cast<afw::coord::CoordSystem>(-1))
1121 {
1122  WcsPersistenceHelper const & keys = WcsPersistenceHelper::get();
1123  if (!record.getSchema().contains(keys.schema)) {
1124  throw LSST_EXCEPT(
1125  afw::table::io::MalformedArchiveError,
1126  "Incorrect schema for Wcs persistence"
1127  );
1128  }
1129  _setWcslibParams();
1130  Eigen::Matrix2d cd = Eigen::Map<Eigen::Matrix2d const>(record[keys.cd].getData());
1131  initWcsLib(
1132  record.get(keys.crval), record.get(keys.crpix), cd,
1133  record.get(keys.ctype1), record.get(keys.ctype2),
1134  record.get(keys.equinox), record.get(keys.radesys),
1135  record.get(keys.cunit1), record.get(keys.cunit2)
1136  );
1137  _initWcs();
1138 }
int _wcsfixCtrl
Do potentially unsafe translations of non-standard unit strings? 0/1 = no/yes.
Definition: Wcs.h:394
table::Key< table::Array< double > > cd
Definition: Wcs.cc:1037
void _setWcslibParams()
Definition: Wcs.cc:69
coord::CoordSystem _coordSystem
Definition: Wcs.h:397
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.
Definition: Wcs.cc:361
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
void _initWcs()
Definition: Wcs.cc:118
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
int _relax
Degree of permissiveness for wcspih (0 for strict); see wcshdr.h for details.
Definition: Wcs.h:393
int _wcshdrCtrl
Controls messages to stderr from wcshdr (0 for none); see wcshdr.h for details.
Definition: Wcs.h:395
Wcs::Wcs ( afwImg::Wcs const &  rhs)
protected

Copy constructor.

Definition at line 454 of file Wcs.cc.

454  :
455  daf::base::Citizen(typeid(this)),
456  _wcsInfo(NULL),
457  _nWcsInfo(rhs._nWcsInfo),
458  _relax(rhs._relax),
459  _wcsfixCtrl(rhs._wcsfixCtrl),
460  _wcshdrCtrl(rhs._wcshdrCtrl),
461  _nReject(rhs._nReject),
462  _coordSystem(static_cast<afwCoord::CoordSystem>(-1))
463 {
464 
465  if (rhs._nWcsInfo > 0) {
466  _wcsInfo = static_cast<struct wcsprm *>(calloc(rhs._nWcsInfo, sizeof(struct wcsprm)));
467  if (_wcsInfo == NULL) {
468  throw LSST_EXCEPT(lsst::pex::exceptions::MemoryError, "Cannot allocate WCS info");
469  }
470 
471  _wcsInfo->flag = -1;
472  int alloc = 1; //Unconditionally allocate memory when calling
473  for (int i = 0; i != rhs._nWcsInfo; ++i) {
474  int status = wcscopy(alloc, &rhs._wcsInfo[i], &_wcsInfo[i]);
475  if (status != 0) {
476  wcsvfree(&i, &_wcsInfo);
477  throw LSST_EXCEPT(lsst::pex::exceptions::MemoryError,
478  (boost::format("Could not copy WCS: wcscopy status = %d : %s") %
479  status % wcs_errmsg[status]).str());
480  }
481  }
482  }
483  _initWcs();
484 }
int _wcsfixCtrl
Do potentially unsafe translations of non-standard unit strings? 0/1 = no/yes.
Definition: Wcs.h:394
coord::CoordSystem _coordSystem
Definition: Wcs.h:397
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
void _initWcs()
Definition: Wcs.cc:118
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
int _relax
Degree of permissiveness for wcspih (0 for strict); see wcshdr.h for details.
Definition: Wcs.h:393
int status
int _wcshdrCtrl
Controls messages to stderr from wcshdr (0 for none); see wcshdr.h for details.
Definition: Wcs.h:395

Member Function Documentation

void Wcs::_initWcs ( )
protected

Definition at line 118 of file Wcs.cc.

119 {
120  if (_wcsInfo) {
122 
123  // tell WCSlib that values have been updated
124  _wcsInfo->flag = 0;
125  // and then tell it to do its internal magic.
126  int status = wcsset(_wcsInfo);
127  if (status != 0) {
128  throw LSST_EXCEPT(except::RuntimeError,
129  (boost::format("Failed to setup wcs structure with wcsset. Status %d: %s") %
130  status % wcs_errmsg[status] ).str());
131  }
132  }
133 }
coord::CoordSystem _coordSystem
Definition: Wcs.h:397
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
int status
CoordSystem makeCoordEnum(std::string const system)
A utility function to get the enum value of a coordinate system from a string name.
Definition: Coord.cc:117
bool Wcs::_isSubset ( Wcs const &  other) const
protectedvirtual

Reimplemented in lsst.afw.image::TanWcs.

Definition at line 520 of file Wcs.cc.

520  {
521  CHECK_NULLS(_wcsInfo, rhs._wcsInfo);
522  CHECK_NULLS(_wcsInfo->crval, rhs._wcsInfo->crval);
523  CHECK_NULLS(_wcsInfo->cd, rhs._wcsInfo->cd);
524  CHECK_NULLS(_wcsInfo->crpix, rhs._wcsInfo->crpix);
525  CHECK_NULLS(_wcsInfo->cunit, rhs._wcsInfo->cunit);
526  CHECK_NULLS(_wcsInfo->ctype, rhs._wcsInfo->ctype);
527  return _nWcsInfo == rhs._nWcsInfo &&
528  _coordSystem == rhs._coordSystem &&
529  _wcsInfo->naxis == rhs._wcsInfo->naxis &&
530  _wcsInfo->equinox == rhs._wcsInfo->equinox &&
531  _wcsInfo->altlin == rhs._wcsInfo->altlin &&
532  compareArrays(_wcsInfo->crval, rhs._wcsInfo->crval, 2) &&
533  compareArrays(_wcsInfo->crpix, rhs._wcsInfo->crpix, 2) &&
534  compareArrays(_wcsInfo->cd, rhs._wcsInfo->cd, 4) &&
535  compareStringArrays(_wcsInfo->cunit, rhs._wcsInfo->cunit, 2) &&
536  compareStringArrays(_wcsInfo->ctype, rhs._wcsInfo->ctype, 2) &&
538  _wcsInfo->crval[1] * afwGeom::degrees) ==
539  rhs.skyToPixel(_wcsInfo->crval[0] * afwGeom::degrees,
540  _wcsInfo->crval[1] * afwGeom::degrees) &&
541  *pixelToSky(_wcsInfo->crpix[0], _wcsInfo->crpix[1]) ==
542  *rhs.pixelToSky(_wcsInfo->crpix[0], _wcsInfo->crpix[1]);
543 }
coord::CoordSystem _coordSystem
Definition: Wcs.h:397
geom::Point2D skyToPixel(geom::Angle sky1, geom::Angle sky2) const
Convert from sky coordinates (e.g. RA/dec) to pixel positions.
Definition: Wcs.cc:782
boost::shared_ptr< coord::Coord > pixelToSky(double pix1, double pix2) const
Convert from pixel position to sky coordinates (e.g. RA/dec)
Definition: Wcs.cc:858
AngleUnit const degrees
Definition: Angle.h:92
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
#define CHECK_NULLS(a, b)
Definition: Wcs.cc:509
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.

1096  {
1097  if (_wcsInfo[0].naxis != 2) return false;
1098  if (std::strcmp(_wcsInfo[0].cunit[0], "deg") != 0) return false;
1099  if (std::strcmp(_wcsInfo[0].cunit[1], "deg") != 0) return false;
1100 
1101  return true;
1102 }
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
void lsst.afw.image::Wcs::_setWcslibParams ( )
protected

Definition at line 69 of file Wcs.cc.

70 {
71  _wcsfixCtrl = // ctrl for wcsfix
72  2; // Translate "H" to "h"
73  _wcshdrCtrl = // ctrl for wcspih
74  2; // Report each rejected keyrecord and the reason why it was rejected
75  _relax = // relax parameter for wcspih;
76  WCSHDR_all; // Accept all extensions recognized by the parser
77 }
int _wcsfixCtrl
Do potentially unsafe translations of non-standard unit strings? 0/1 = no/yes.
Definition: Wcs.h:394
int _relax
Degree of permissiveness for wcspih (0 for strict); see wcshdr.h for details.
Definition: Wcs.h:393
int _wcshdrCtrl
Controls messages to stderr from wcshdr (0 for none); see wcshdr.h for details.
Definition: Wcs.h:395
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.

552  {
553  return Wcs::Ptr(new Wcs(*this));
554 }
boost::shared_ptr< Wcs > Ptr
Definition: Wcs.h:113
Wcs()
Construct an invalid Wcs given no arguments.
Definition: Wcs.cc:88
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.

778  {
779  return coord.convert(_coordSystem);
780 }
coord::CoordSystem _coordSystem
Definition: Wcs.h:397
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.

591  {
592  assert(_wcsInfo);
593 
594  int const naxis = _wcsInfo->naxis;
595 
596  //If naxis != 2, I'm not sure if any of what follows is correct
597  assert(naxis == 2);
598  //Origin is at (1,1). Adjust to avoid off by one.
599  if (flipLR) {
600  _wcsInfo->cd[0] = -_wcsInfo->cd[0];
601  _wcsInfo->cd[2] = -_wcsInfo->cd[2];
602  _wcsInfo->crpix[0] = -_wcsInfo->crpix[0] + dimensions.getX() + 1;
603  }
604  if (flipTB) {
605  _wcsInfo->cd[1] = -_wcsInfo->cd[1];
606  _wcsInfo->cd[3] = -_wcsInfo->cd[3];
607  _wcsInfo->crpix[1] = -_wcsInfo->crpix[1]+dimensions.getY() + 1;
608  }
609 
610  // tells libwcs to invalidate cached data, since transformation has been modified
611  _wcsInfo->flag = 0;
612 }
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
Eigen::Matrix2d Wcs::getCDMatrix ( ) const

Returns the CD matrix.

Definition at line 573 of file Wcs.cc.

573  {
574  assert(_wcsInfo);
575  int const naxis = _wcsInfo->naxis;
576 
577  //If naxis != 2, I'm not sure if any of what follows is correct
578  assert(naxis == 2);
579 
580  Eigen::Matrix2d C;
581 
582  for (int i=0; i< naxis; ++i){
583  for (int j=0; j<naxis; ++j) {
584  C(i,j) = _wcsInfo->cd[ (i*naxis) + j ];
585  }
586  }
587 
588  return C;
589 }
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
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.

668  {
670 }
static lsst::daf::base::PropertyList::Ptr generatePropertySet(lsst::afw::image::Wcs const &wcs)
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.

1009 {
1011 }
A 2D linear coordinate transformation.
Eigen::Matrix2d getCDMatrix() const
Returns the CD matrix.
Definition: Wcs.cc:573
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.

1074 { return getWcsPersistenceName(); }
GeomPoint Wcs::getPixelOrigin ( ) const

Returns CRPIX (corrected to LSST convention).

Definition at line 565 of file Wcs.cc.

565  {
566  assert(_wcsInfo);
567  //Convert from fits units back to lsst units
568  double p1 = _wcsInfo->crpix[0] + fitsToLsstPixels;
569  double p2 = _wcsInfo->crpix[1] + fitsToLsstPixels;
570  return afwGeom::Point2D(p1, p2);
571 }
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
const int fitsToLsstPixels
Definition: Wcs.cc:80
Point< double, 2 > Point2D
Definition: Point.h:286
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.

1076 { 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.

560  {
561  assert(_wcsInfo);
562  return makeCorrectCoord(_wcsInfo->crval[0] * afwGeom::degrees, _wcsInfo->crval[1] * afwGeom::degrees);
563 }
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 C...
Definition: Wcs.cc:879
AngleUnit const degrees
Definition: Angle.h:92
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
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.

220 { 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
crvalThe sky position of the reference point
crpixThe pixel position corresponding to crval in LSST units
CDMatrix describing transformations from pixel to sky positions
ctype1Projection system used (see description of Wcs)
ctype2Projection system used (see description of Wcs)
equinoxEquinox of coordinate system, eg 2000 (Julian) or 1950 (Besselian)
raDecSysSystem used to describe right ascension or declination, e.g FK4, FK5 or ICRS
cunits1Units of sky position. One of deg, arcmin or arcsec
cunits2Units of sky position. One of deg, arcmin or arcsec

Definition at line 361 of file Wcs.cc.

364  {
365 
366  //Check CD is a valid size
367  if( (CD.rows() != 2) || (CD.cols() != 2) ) {
368  string msg = "CD is not a 2x2 matrix";
369  throw LSST_EXCEPT(except::InvalidParameterError, msg);
370  }
371 
372  //Check that cunits are legitimate values
373  bool isValid = (cunits1 == "deg");
374  isValid |= (cunits1 == "arcmin");
375  isValid |= (cunits1 == "arcsec");
376  isValid |= (cunits1 == "mas");
377 
378  if (!isValid) {
379  string msg = "CUNITS1 must be one of {deg|arcmin|arcsec|mas}";
380  throw LSST_EXCEPT(except::InvalidParameterError, msg);
381  }
382 
383  isValid = (cunits2 == "deg");
384  isValid |= (cunits2 == "arcmin");
385  isValid |= (cunits2 == "arcsec");
386  isValid |= (cunits2 == "mas");
387 
388  if (!isValid) {
389  string msg = "CUNITS2 must be one of {deg|arcmin|arcsec|mas}";
390  throw LSST_EXCEPT(except::InvalidParameterError, msg);
391  }
392 
393  //Initialise the wcs struct
394  _wcsInfo = static_cast<struct wcsprm *>(malloc(sizeof(struct wcsprm)));
395  if (_wcsInfo == NULL) {
396  throw LSST_EXCEPT(except::MemoryError, "Cannot allocate WCS info");
397  }
398 
399  _wcsInfo->flag = -1;
400  int status = wcsini(true, 2, _wcsInfo); //2 indicates a naxis==2, a two dimensional image
401  if(status != 0) {
402  throw LSST_EXCEPT(except::MemoryError,
403  (boost::format("Failed to allocate memory with wcsini. Status %d: %s") %
404  status % wcs_errmsg[status] ).str());
405  }
406 
407  //Set crval, crpix and CD. Internally to the class, we use fits units for consistency with
408  //wcslib.
409  _wcsInfo->crval[0] = crval.getX();
410  _wcsInfo->crval[1] = crval.getY();
411  _wcsInfo->crpix[0] = crpix.getX() + lsstToFitsPixels;
412  _wcsInfo->crpix[1] = crpix.getY() + lsstToFitsPixels;
413 
414  //Set the CD matrix
415  for (int i=0; i<2; ++i) {
416  for (int j=0; j<2; ++j) {
417  _wcsInfo->cd[(2*i) + j] = CD(i,j);
418  }
419  }
420 
421  //Specify that we have a CD matrix, but no PC or CROTA
422  _wcsInfo->altlin = 2;
423  _wcsInfo->flag = 0; //values have been updated
424 
425  //This is a work around for what I think is a bug in wcslib. ->types is neither
426  //initialised or set to NULL by default, so if I try to delete a Wcs object,
427  //wcslib then attempts to free non-existent space, and the code can crash.
428  _wcsInfo->types = NULL;
429 
430  //Set the coordinate system
431  strncpy(_wcsInfo->ctype[0], ctype1.c_str(), STRLEN);
432  strncpy(_wcsInfo->ctype[1], ctype2.c_str(), STRLEN);
433  strncpy(_wcsInfo->radesys, raDecSys.c_str(), STRLEN);
434  _wcsInfo->equinox = equinox;
435 
436  //Set the units
437  strncpy(_wcsInfo->cunit[0], cunits1.c_str(), STRLEN);
438  strncpy(_wcsInfo->cunit[1], cunits2.c_str(), STRLEN);
439 
440  _nWcsInfo = 1; //Specify that we have only one coordinate representation
441 
442  //Tell wcslib that we are need to set up internal values
443  status=wcsset(_wcsInfo);
444  if(status != 0) {
445  throw LSST_EXCEPT(except::RuntimeError,
446  (boost::format("Failed to setup wcs structure with wcsset. Status %d: %s") %
447  status % wcs_errmsg[status] ).str());
448 
449  }
450 }
const int lsstToFitsPixels
Definition: Wcs.cc:79
table::PointKey< double > crval
Definition: Wcs.cc:1035
table::Key< std::string > ctype2
Definition: Wcs.cc:1039
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
table::Key< double > equinox
Definition: Wcs.cc:1040
const int STRLEN
Definition: Wcs.cc:66
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
int status
table::Key< std::string > ctype1
Definition: Wcs.cc:1038
table::PointKey< double > crpix
Definition: Wcs.cc:1036
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.

173  {
178  class HeaderAccess {
179  public:
181  CONST_PTR(lsst::daf::base::PropertySet) const& toRead() { return _constHeader; }
183  PTR(lsst::daf::base::PropertySet) const& toWrite() {
184  if (!_hackHeader) {
185  _hackHeader = _constHeader->deepCopy();
186  _constHeader = _hackHeader;
187  }
188  return _hackHeader;
189  }
190 
192  HeaderAccess(CONST_PTR(lsst::daf::base::PropertySet) const& header) :
193  _constHeader(header), _hackHeader() {}
194 
195  private:
197  PTR(lsst::daf::base::PropertySet) _hackHeader;
198  };
199 
200  HeaderAccess access(header);
201 
202  // Some headers (e.g. SDSS ones from FNAL) have EQUINOX as a string. Fix this,
203  // as wcslib 4.4.4 refuses to handle it. Furthermore, some headers (e.g. Skymapper)
204  // use a string but prepend 'J': "J2000.0" -- if we don't handle this, we deduce
205  // an equinox of 0.0
206  {
207  std::string const& key = "EQUINOX";
208  if (access.toRead()->exists(key) && access.toRead()->typeOf(key) == typeid(std::string)) {
209  auto equinoxStr = access.toRead()->getAsString(key).c_str();
210  double equinox = ::atof(equinoxStr[0] == 'J' ? equinoxStr + 1 : equinoxStr);
211  access.toWrite()->set(key, equinox);
212  }
213  }
214 
215  //Check header isn't empty
216  int nCards = lsst::afw::formatters::countFitsHeaderCards(access.toRead());
217  if (nCards <= 0) {
218  string msg = "Could not parse FITS WCS: no header cards found";
219  throw LSST_EXCEPT(except::InvalidParameterError, msg);
220  }
221 
222  //While the standard does not insist on CRVAL and CRPIX being present, it
223  //is almost certain their absence indicates a problem.
224  //Check for CRPIX
225  if( !access.toRead()->exists("CRPIX1") && !access.toRead()->exists("CRPIX1a")) {
226  string msg = "Neither CRPIX1 not CRPIX1a found";
227  throw LSST_EXCEPT(except::InvalidParameterError, msg);
228  }
229 
230  if( !access.toRead()->exists("CRPIX2") && !access.toRead()->exists("CRPIX2a")) {
231  string msg = "Neither CRPIX2 not CRPIX2a found";
232  throw LSST_EXCEPT(except::InvalidParameterError, msg);
233  }
234 
235  //And the same for CRVAL
236  if( !access.toRead()->exists("CRVAL1") && !access.toRead()->exists("CRVAL1a")) {
237  string msg = "Neither CRVAL1 not CRVAL1a found";
238  throw LSST_EXCEPT(except::InvalidParameterError, msg);
239  }
240 
241  if( !access.toRead()->exists("CRVAL2") && !access.toRead()->exists("CRVAL2a")) {
242  string msg = "Neither CRVAL2 not CRVAL2a found";
243  throw LSST_EXCEPT(except::InvalidParameterError, msg);
244  }
245  /*
246  * According to Greisen and Calabretta (A&A 395, 1061–1075 (2002)) it's illegal to mix PCi_j and CDi_j
247  * headers; unfortunately Subaru puts both in its headers. It actually uses PC001002 instead of PC1_2
248  * (dating to a proposed FITS standard from 1996) and at least sometimes fails to include CDELT[12],
249  * so the CD and PC matrices are inconsistent
250  *
251  * If we detect any part of a CD matrix, delete all PC matrices
252  */
253  if(access.toRead()->exists("CD1_1") || access.toRead()->exists("CD1_2") ||
254  access.toRead()->exists("CD2_1") || access.toRead()->exists("CD2_2")) {
255  for (int i = 1; i <= 2; ++i) {
256  for (int j = 1; j <= 2; ++j) {
257  std::string key = (boost::format("PC%i_%i") % j % i).str();
258  if (access.toRead()->exists(key)) {
259  double const val = access.toRead()->getAsDouble(key);
260  access.toWrite()->remove(key);
261  access.toWrite()->add("X_" + key, val);
262  }
263 
264  key = (boost::format("PC%03d%03d") % j % i).str();
265  if (access.toRead()->exists(key)) {
266  double const val = access.toRead()->getAsDouble(key);
267  access.toWrite()->remove(key);
268  access.toWrite()->add("X_" + key, val);
269  }
270  }
271  }
272  }
273 
274  //Pass the header into wcslib's formatter to extract & setup the Wcs. First need
275  //to convert to a C style string, so the compile doesn't complain about constness
276  std::string metadataStr = lsst::afw::formatters::formatFitsProperties(access.toRead());
277  // We own the data, and wcslib is slack about constness, so no qualms with casting away const
278  char *hdrString = const_cast<char*>(metadataStr.c_str());
279  //printf("wcspih string:\n%s\n", hdrString);
280 
281  nCards = lsst::afw::formatters::countFitsHeaderCards(access.toRead()); // we may have dropped some
282  int pihStatus = wcspih(hdrString, nCards, _relax, _wcshdrCtrl, &_nReject, &_nWcsInfo, &_wcsInfo);
283 
284  if (pihStatus != 0) {
285  throw LSST_EXCEPT(except::RuntimeError,
286  (boost::format("Could not parse FITS WCS: wcspih status = %d (%s)") %
287  pihStatus % wcs_errmsg[pihStatus]).str());
288  }
289 
290  //Run wcsfix on _wcsInfo to try and fix any problems it knows about.
291  const int *naxes = NULL; // should be {NAXIS1, NAXIS2, ...} to check cylindrical projections
292  int stats[NWCSFIX]; // status returns from wcsfix
293  int fixStatus = wcsfix(_wcsfixCtrl, naxes, _wcsInfo, stats);
294  if (fixStatus != 0) {
295  std::stringstream errStream;
296  errStream << "Could not parse FITS WCS: wcsfix failed " << std::endl;
297  for (int ii = 0; ii < NWCSFIX; ++ii) {
298  if (stats[ii] >= 0) {
299  errStream << "\t" << ii << ": " << stats[ii] << " " << wcsfix_errmsg[stats[ii]] << std::endl;
300  } else {
301  errStream << "\t" << ii << ": " << stats[ii] << std::endl;
302  }
303  }
304  }
305 
306  //The Wcs standard requires a default value for RADESYS if the keyword
307  //doesn't exist in header, but wcslib doesn't set it. So we do so here. This code
308  //conforms to Calabretta & Greisen 2002 \S 3.1
309  if (!(access.toRead()->exists("RADESYS") || access.toRead()->exists("RADESYSa"))) {
310  // If RADECSYS exists, use that (counter to Calabretta & Greisen 2002 \S 3.1, but commonly used).
311  // If equinox exist and < 1984, use FK4. If >= 1984, use FK5
312  if (access.toRead()->exists("RADECSYS")) {
313  strncpy(_wcsInfo->radesys, access.toRead()->getAsString("RADECSYS").c_str(), STRLEN);
314  } else if (access.toRead()->exists("EQUINOX") || access.toRead()->exists("EQUINOXa")) {
315  std::string const EQUINOX = access.toRead()->exists("EQUINOX") ? "EQUINOX" : "EQUINOXa";
316  double const equinox = access.toRead()->getAsDouble(EQUINOX);
317  if(equinox < 1984) {
318  strncpy(_wcsInfo->radesys, "FK4", STRLEN);
319  } else {
320  strncpy(_wcsInfo->radesys, "FK5", STRLEN);
321  }
322  } else {
323  //If Equinox doesn't exist, default to ICRS
324  strncpy(_wcsInfo->radesys, "ICRS", STRLEN);
325  }
326  }
327  // strip trailing whitespace
328  {
329  for(int i = strlen(_wcsInfo->radesys) - 1; i >= 0; i--) {
330  if (isspace(_wcsInfo->radesys[i])) {
331  _wcsInfo->radesys[i] = '\0';
332  }
333  }
334  }
335  //
336  // If there are no CDi_j cards in the header, set CDi_j from PCi_j
337  // CDi_j == CDELTi*PCi_j
338  //
339  if ((_wcsInfo->altlin & 2) == 0) { // no CDi_j cards were found in the header
340  double const *cdelt = _wcsInfo->cdelt;
341  double const *pc = _wcsInfo->pc;
342  double *cd = _wcsInfo->cd;
343 
344  cd[0] = cdelt[0]*pc[0]; // 1_1
345  cd[1] = cdelt[0]*pc[1]; // 1_2
346  cd[2] = cdelt[1]*pc[2]; // 2_1
347  cd[3] = cdelt[1]*pc[3]; // 2_2
348  }
349 }
int _wcsfixCtrl
Do potentially unsafe translations of non-standard unit strings? 0/1 = no/yes.
Definition: Wcs.h:394
#define PTR(...)
Definition: base.h:41
table::Key< table::Array< double > > cd
Definition: Wcs.cc:1037
#define CONST_PTR(...)
Definition: base.h:47
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
table::Key< double > equinox
Definition: Wcs.cc:1040
const int STRLEN
Definition: Wcs.cc:66
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
int _relax
Degree of permissiveness for wcspih (0 for strict); see wcshdr.h for details.
Definition: Wcs.h:393
Class for storing generic metadata.
Definition: PropertySet.h:82
std::string formatFitsProperties(boost::shared_ptr< lsst::daf::base::PropertySet const > const &prop)
int _wcshdrCtrl
Controls messages to stderr from wcshdr (0 for none); see wcshdr.h for details.
Definition: Wcs.h:395
ImageT val
Definition: CR.cc:154
int countFitsHeaderCards(boost::shared_ptr< lsst::daf::base::PropertySet const > const &prop)
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.

672  {
673  // We calculate the determinant of the CD (i.e the rotation and scaling)
674  // matrix. If this determinant is positive, then the image can be rotated
675  // to a position where increasing the right ascension and declination
676  // increases the horizontal and vertical pixel position. In this case the
677  // image is flipped.
678  assert(_wcsInfo);
679  double det = (_wcsInfo->cd[0] * _wcsInfo->cd[3]) - (_wcsInfo->cd[1] * _wcsInfo->cd[2]);
680 
681  if (det == 0) {
682  throw(LSST_EXCEPT(except::RuntimeError, "Wcs CD matrix is singular"));
683  }
684 
685  return (det > 0);
686 }
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
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.

1104  {
1105  if (!_mayBePersistable()) {
1106  return false;
1107  }
1108  // The current table persistence only works for TAN and TAN-SIP projections
1109  return false;
1110 }
bool _mayBePersistable() const
Perform basic checks on whether *this might be persistable.
Definition: Wcs.cc:1096
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]coordPosition in sky coordinates where transform is desired.
[in]skyUnitUnits to use for sky coordinates; units of matrix elements will be skyUnits/pixel.

Definition at line 934 of file Wcs.cc.

937  {
938  return linearizePixelToSkyInternal(skyToPixel(coord), coord, skyUnit);
939 }
geom::Point2D skyToPixel(geom::Angle sky1, geom::Angle sky2) const
Convert from sky coordinates (e.g. RA/dec) to pixel positions.
Definition: Wcs.cc:782
virtual geom::AffineTransform linearizePixelToSkyInternal(geom::Point2D const &pix, coord::Coord const &coord, geom::AngleUnit skyUnit) const
Definition: Wcs.cc:951
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]pixPosition in pixel coordinates where transform is desired.
[in]skyUnitUnits to use for sky coordinates; units of matrix elements will be skyUnits/pixel.

Definition at line 940 of file Wcs.cc.

943  {
944  return linearizePixelToSkyInternal(pix, *pixelToSky(pix), skyUnit);
945 }
boost::shared_ptr< coord::Coord > pixelToSky(double pix1, double pix2) const
Convert from pixel position to sky coordinates (e.g. RA/dec)
Definition: Wcs.cc:858
virtual geom::AffineTransform linearizePixelToSkyInternal(geom::Point2D const &pix, coord::Coord const &coord, geom::AngleUnit skyUnit) const
Definition: Wcs.cc:951
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.

955  {
956  //
957  // Figure out the (0, 0), (0, 1), and (1, 0) ra/dec coordinates of the corners of a square drawn in pixel
958  // It'd be better to centre the square at sky00, but that would involve another conversion between sky and
959  // pixel coordinates so I didn't bother
960  //
961  const double side = 10; // length of the square's sides in pixels
962  GeomPoint const sky00 = coord.getPosition(skyUnit);
963  typedef std::pair<lsst::afw::geom::Angle, lsst::afw::geom::Angle> AngleAngle;
964  AngleAngle const dsky10 = coord.getTangentPlaneOffset(*pixelToSky(pix00 + afwGeom::Extent2D(side, 0)));
965  AngleAngle const dsky01 = coord.getTangentPlaneOffset(*pixelToSky(pix00 + afwGeom::Extent2D(0, side)));
966 
967  Eigen::Matrix2d m;
968  m(0, 0) = dsky10.first.asAngularUnits(skyUnit)/side;
969  m(0, 1) = dsky01.first.asAngularUnits(skyUnit)/side;
970  m(1, 0) = dsky10.second.asAngularUnits(skyUnit)/side;
971  m(1, 1) = dsky01.second.asAngularUnits(skyUnit)/side;
972 
973  Eigen::Vector2d sky00v;
974  sky00v << sky00.getX(), sky00.getY();
975  Eigen::Vector2d pix00v;
976  pix00v << pix00.getX(), pix00.getY();
977  //return lsst::afw::geom::AffineTransform(m, lsst::afw::geom::Extent2D(sky00v - m * pix00v));
978  return lsst::afw::geom::AffineTransform(m, (sky00v - m * pix00v));
979 }
boost::shared_ptr< coord::Coord > pixelToSky(double pix1, double pix2) const
Convert from pixel position to sky coordinates (e.g. RA/dec)
Definition: Wcs.cc:858
An affine coordinate transformation consisting of a linear transformation and an offset.
tuple m
Definition: lsstimport.py:48
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]coordPosition in sky coordinates where transform is desired.
[in]skyUnitUnits to use for sky coordinates; units of matrix elements will be pixels/skyUnit.

Definition at line 981 of file Wcs.cc.

984  {
985  return linearizeSkyToPixelInternal(skyToPixel(coord), coord, skyUnit);
986 }
geom::Point2D skyToPixel(geom::Angle sky1, geom::Angle sky2) const
Convert from sky coordinates (e.g. RA/dec) to pixel positions.
Definition: Wcs.cc:782
virtual geom::AffineTransform linearizeSkyToPixelInternal(geom::Point2D const &pix, coord::Coord const &coord, geom::AngleUnit skyUnit) const
Definition: Wcs.cc:999
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]pixPosition in pixel coordinates where transform is desired.
[in]skyUnitUnits to use for sky coordinates; units of matrix elements will be pixels/skyUnit.

Definition at line 988 of file Wcs.cc.

991  {
992  return linearizeSkyToPixelInternal(pix, *pixelToSky(pix), skyUnit);
993 }
virtual geom::AffineTransform linearizeSkyToPixelInternal(geom::Point2D const &pix, coord::Coord const &coord, geom::AngleUnit skyUnit) const
Definition: Wcs.cc:999
boost::shared_ptr< coord::Coord > pixelToSky(double pix1, double pix2) const
Convert from pixel position to sky coordinates (e.g. RA/dec)
Definition: Wcs.cc:858
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.

1003  {
1004  lsst::afw::geom::AffineTransform inverse = linearizePixelToSkyInternal(pix00, coord, skyUnit);
1005  return inverse.invert();
1006 }
AffineTransform const invert() const
An affine coordinate transformation consisting of a linear transformation and an offset.
virtual geom::AffineTransform linearizePixelToSkyInternal(geom::Point2D const &pix, coord::Coord const &coord, geom::AngleUnit skyUnit) const
Definition: Wcs.cc:951
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.

879  {
880 
881  //Construct a coord object of the correct type
882  int const ncompare = 4; // we only care about type's first 4 chars
883  char *type = _wcsInfo->ctype[0];
884  char *radesys = _wcsInfo->radesys;
885  double equinox = _wcsInfo->equinox;
886 
887  if (strncmp(type, "RA--", ncompare) == 0) { // Our default. If it's often something else, consider
888  ; // using an tr1::unordered_map
889  if(strcmp(radesys, "ICRS") == 0) {
890  return afwCoord::makeCoord(afwCoord::ICRS, sky0, sky1);
891  }
892  if(strcmp(radesys, "FK5") == 0) {
893  return afwCoord::makeCoord(afwCoord::FK5, sky0, sky1, equinox);
894  } else {
895  throw LSST_EXCEPT(except::RuntimeError,
896  (boost::format("Can't create Coord object: Unrecognised radesys %s") %
897  radesys).str());
898  }
899 
900  } else if (strncmp(type, "GLON", ncompare) == 0) {
901  return afwCoord::makeCoord(afwCoord::GALACTIC, sky0, sky1);
902  } else if (strncmp(type, "ELON", ncompare) == 0) {
903  return afwCoord::makeCoord(afwCoord::ECLIPTIC, sky0, sky1, equinox);
904  } else if (strncmp(type, "DEC-", ncompare) == 0) {
905  //check for the case where the ctypes are swapped. Note how sky0 and sky1 are swapped as well
906 
907  //Our default
908  if(strcmp(radesys, "ICRS") == 0) {
909  return afwCoord::makeCoord(afwCoord::ICRS, sky1, sky0);
910  }
911  if(strcmp(radesys, "FK5") == 0) {
912  return afwCoord::makeCoord(afwCoord::FK5, sky1, sky0, equinox);
913  } else {
914  throw LSST_EXCEPT(except::RuntimeError,
915  (boost::format("Can't create Coord object: Unrecognised radesys %s") %
916  radesys).str());
917  }
918  } else if (strncmp(type, "GLAT", ncompare) == 0) {
919  return afwCoord::makeCoord(afwCoord::GALACTIC, sky1, sky0);
920  } else if (strncmp(type, "ELAT", ncompare) == 0) {
921  return afwCoord::makeCoord(afwCoord::ECLIPTIC, sky1, sky0, equinox);
922  } else {
923  //Give up in disgust
924  throw LSST_EXCEPT(except::RuntimeError,
925  (boost::format("Can't create Coord object: Unrecognised sys %s") %
926  type).str());
927  }
928 
929  //Can't get here
930  assert(0);
931 }
table::Key< std::string > radesys
Definition: Wcs.cc:1041
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
table::Key< double > equinox
Definition: Wcs.cc:1040
Coord::Ptr makeCoord(CoordSystem const system, lsst::afw::geom::Angle const ra, lsst::afw::geom::Angle const dec, double const epoch)
Factory function to create a Coord of arbitrary type with decimal RA,Dec.
Definition: Coord.cc:1211
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
bool lsst.afw.image::Wcs::operator!= ( Wcs const &  other) const
inline

Definition at line 135 of file Wcs.h.

135 { 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.

486  {
487  if (&other == this) return true;
488  // We do a bidirectional test with a virtual member function in case one of us is a derived
489  // class with members we don't know about here.
490  // This is not the most efficient possible implementation, but I think it's the easiest one
491  // with which to ensure correctness, and I think that's more important in this case.
492  return this->_isSubset(other) && other._isSubset(*this);
493 }
virtual bool _isSubset(Wcs const &other) const
Definition: Wcs.cc:520
double Wcs::pixArea ( lsst::afw::geom::Point2D  pix00) const

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

Parameters
pix00The pixel point where the area is desired

Definition at line 692 of file Wcs.cc.

693  {
694  //
695  // Figure out the (0, 0), (0, 1), and (1, 0) ra/dec coordinates of the corners of a square drawn in pixel
696  // It'd be better to centre the square at sky00, but that would involve another conversion between sky and
697  // pixel coordinates so I didn't bother
698  //
699  const double side = 1; // length of the square's sides in pixels
700 
701  // Work in 3-space to avoid RA wrapping and pole issues.
702  afwGeom::Point3D v0 = pixelToSky(pix00)->getVector();
703 
704  // Step by "side" in x and y pixel directions...
705  GeomPoint px(pix00);
706  GeomPoint py(pix00);
707  px.shift(afwGeom::Extent2D(side, 0));
708  py.shift(afwGeom::Extent2D(0, side));
709  // Push the points through the WCS, and find difference in 3-space.
710  afwGeom::Extent3D dx = pixelToSky(px)->getVector() - v0;
711  afwGeom::Extent3D dy = pixelToSky(py)->getVector() - v0;
712 
713  // Compute |cross product| = area of parallelogram with sides dx,dy
714  // FIXME -- this is slightly incorrect -- it's making the small-angle
715  // approximation, taking the distance *through* the unit sphere
716  // rather than over its surface.
717  // This is in units of ~radians^2
718  double area = sqrt(square(dx[1]*dy[2] - dx[2]*dy[1]) +
719  square(dx[2]*dy[0] - dx[0]*dy[2]) +
720  square(dx[0]*dy[1] - dx[1]*dy[0]));
721 
722  return area / square(side) * square(180. / afwGeom::PI);
723 }
boost::shared_ptr< coord::Coord > pixelToSky(double pix1, double pix2) const
Convert from pixel position to sky coordinates (e.g. RA/dec)
Definition: Wcs.cc:858
A coordinate class intended to represent absolute positions.
Definition: PSF.h:39
double const PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:18
afwGeom::Angle Wcs::pixelScale ( ) const

Returns the pixel scale [Angle/pixel].

Definition at line 725 of file Wcs.cc.

725  {
726  return sqrt(pixArea(getPixelOrigin())) * afwGeom::degrees;
727 }
AngleUnit const degrees
Definition: Angle.h:92
lsst::afw::geom::Point2D getPixelOrigin() const
Returns CRPIX (corrected to LSST convention).
Definition: Wcs.cc:565
double pixArea(lsst::afw::geom::Point2D pix00) const
Sky area covered by a pixel at position pix00 in units of square degrees.
Definition: Wcs.cc:692
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.

858  {
859  assert(_wcsInfo);
860 
861  afwGeom::Angle skyTmp[2];
862  pixelToSkyImpl(pixel1, pixel2, skyTmp);
863  return makeCorrectCoord(skyTmp[0], skyTmp[1]);
864 }
virtual void pixelToSkyImpl(double pixel1, double pixel2, geom::Angle skyTmp[2]) const
Definition: Wcs.cc:831
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 C...
Definition: Wcs.cc:879
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
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.

854  {
855  return pixelToSky(pixel.getX(), pixel.getY());
856 }
boost::shared_ptr< coord::Coord > pixelToSky(double pix1, double pix2) const
Convert from pixel position to sky coordinates (e.g. RA/dec)
Definition: Wcs.cc:858
table::PointKey< int > pixel
void Wcs::pixelToSky ( double  pixel1,
double  pixel2,
geom::Angle sky1,
geom::Angle sky2 
) const

Definition at line 866 of file Wcs.cc.

866  {
867  afwGeom::Angle skyTmp[2];
868  // HACK -- we shouldn't need to initialize these -- pixelToSkyImpl() sets them unless an
869  // exception is thrown -- but be safe.
870  skyTmp[0] = 0. * afwGeom::radians;
871  skyTmp[1] = 0. * afwGeom::radians;
872  pixelToSkyImpl(pixel1, pixel2, skyTmp);
873  sky1 = skyTmp[0];
874  sky2 = skyTmp[1];
875 }
AngleUnit const radians
constant with units of radians
Definition: Angle.h:91
virtual void pixelToSkyImpl(double pixel1, double pixel2, geom::Angle skyTmp[2]) const
Definition: Wcs.cc:831
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.

832 {
833  assert(_wcsInfo);
834 
835  // wcslib assumes 1-indexed coordinates
836  double pixTmp[2] = { pixel1 - lsst::afw::image::PixelZeroPos + lsstToFitsPixels,
837  pixel2 - lsst::afw::image::PixelZeroPos + lsstToFitsPixels};
838  double imgcrd[2];
839  double phi, theta;
840 
841  double sky[2];
842  int status = 0;
843  status = wcsp2s(_wcsInfo, 1, 2, pixTmp, imgcrd, &phi, &theta, sky, &status);
844  if (status > 0) {
845  throw LSST_EXCEPT(except::RuntimeError,
846  (boost::format("Error: wcslib returned a status code of %d at pixel %s, %s: %s") %
847  status % pixel1 % pixel2 % wcs_errmsg[status]).str());
848  }
849  // FIXME -- _wcsInfo.lat, _wcsInfo.lng ?
850  skyTmp[0] = sky[0] * afwGeom::degrees;
851  skyTmp[1] = sky[1] * afwGeom::degrees;
852 }
const int lsstToFitsPixels
Definition: Wcs.cc:79
AngleUnit const degrees
Definition: Angle.h:92
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
int status
const double PixelZeroPos
FITS uses 1.0, SDSS uses 0.5, LSST uses 0.0 (http://dev.lsstcorp.org/trac/wiki/BottomLeftPixelProposa...
Definition: ImageUtils.h:42
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.

614  {
615  assert(_wcsInfo);
616 
617  while (nQuarter < 0 ) {
618  nQuarter += 4;
619  }
620 
621 
622  int const naxis = _wcsInfo->naxis;
623 
624  //If naxis != 2, I'm not sure if any of what follows is correct
625  assert(naxis == 2);
626  double a = _wcsInfo->cd[0];
627  double b = _wcsInfo->cd[1];
628  double c = _wcsInfo->cd[2];
629  double d = _wcsInfo->cd[3];
630  double crpx = _wcsInfo->crpix[0];
631  double crpy = _wcsInfo->crpix[1];
632  //Origin is at (1,1). Adjust to avoid off by one.
633  //E.g. CRPIX one pixel off the UR of the grid should go to
634  //(0,0) for nQuarter=2
635  switch (nQuarter%4) {
636  case 0:
637  break;
638  case 1:
639  _wcsInfo->cd[0] = -b;
640  _wcsInfo->cd[1] = a;
641  _wcsInfo->cd[2] = -d;
642  _wcsInfo->cd[3] = c;
643  _wcsInfo->crpix[0] = -crpy + dimensions.getY() + 1;
644  _wcsInfo->crpix[1] = crpx;
645  break;
646  case 2:
647  _wcsInfo->cd[0] = -a;
648  _wcsInfo->cd[1] = -b;
649  _wcsInfo->cd[2] = -c;
650  _wcsInfo->cd[3] = -d;
651  _wcsInfo->crpix[0] = -crpx + dimensions.getX() + 1;
652  _wcsInfo->crpix[1] = -crpy + dimensions.getY() + 1;
653  break;
654  case 3:
655  _wcsInfo->cd[0] = b;
656  _wcsInfo->cd[1] = -a;
657  _wcsInfo->cd[2] = d;
658  _wcsInfo->cd[3] = -c;
659  _wcsInfo->crpix[0] = crpy;
660  _wcsInfo->crpix[1] = -crpx + dimensions.getX() + 1;
661  break;
662  }
663 
664  // tells libwcs to invalidate cached data, since transformation has been modified
665  _wcsInfo->flag = 0;
666 }
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
afw::table::Key< double > b
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.

1164  {
1165  assert(_wcsInfo);
1166  _wcsInfo->crpix[0] += dx;
1167  _wcsInfo->crpix[1] += dy;
1168 
1169  // tells libwcs to invalidate cached data, since transformation has been modified
1170  _wcsInfo->flag = 0;
1171 }
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
virtual void lsst.afw.image::Wcs::shiftReferencePixel ( geom::Extent2D const &  d)
inlinevirtual

Definition at line 321 of file Wcs.h.

321 { shiftReferencePixel(d.getX(), d.getY()); }
virtual void shiftReferencePixel(double dx, double dy)
Move the pixel reference position by (dx, dy)
Definition: Wcs.cc:1164
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.

786  {
787  assert(_wcsInfo);
788 
790  double skyTmp[2];
791  double imgcrd[2];
792  double phi, theta;
793  double pixTmp[2];
794 
795  /*
796  printf("skyToIWC: _coordSystem = %i\n", (int)_coordSystem);
797  printf("coord (%.3f, %.3f)\n", coord->getLongitude().asDegrees(), coord->getLatitude().asDegrees());
798  printf("->sky (%.3f, %.3f)\n", sky->getLongitude().asDegrees(), sky->getLatitude().asDegrees());
799  */
800 
801  skyTmp[_wcsInfo->lng] = sky->getLongitude().asDegrees();
802  skyTmp[_wcsInfo->lat] = sky->getLatitude() .asDegrees();
803 
804  //Estimate pixel coordinates
805  int stat[1];
806  int status = 0;
807  imgcrd[0] = imgcrd[1] = -1e6;
808  /*
809  printf(" skyTmp[] = (%.3f, %.3f)\n", skyTmp[0], skyTmp[1]);
810  printf(" _wcsInfo->lng,lat = %i, %i\n", _wcsInfo->lng, _wcsInfo->lat);
811  */
812  status = wcss2p(_wcsInfo, 1, 2, skyTmp, &phi, &theta, imgcrd, pixTmp, stat);
813  if (status > 0) {
814  throw LSST_EXCEPT(except::RuntimeError,
815  (boost::format("Error: wcslib returned a status code of %d at sky %s, %s deg: %s") %
816  status % skyTmp[0] % skyTmp[1] % wcs_errmsg[status]).str());
817  }
818  /*
819  printf("->iwc (%.3f, %.3f)\n", imgcrd[0], imgcrd[1]);
820  printf("-> pix (%.2f, %.2f)\n", pixTmp[0], pixTmp[1]);
821  afwCoord::Coord::Ptr crval = getSkyOrigin();
822  printf("(crval is (%.3f, %.3f))\n", crval->getLongitude().asDegrees(), crval->getLatitude().asDegrees());
823  */
824  return GeomPoint(imgcrd[0], imgcrd[1]);
825 }
boost::shared_ptr< Coord > Ptr
Definition: Coord.h:72
afw::coord::Coord::Ptr convertCoordToSky(coord::Coord const &coord) const
Definition: Wcs.cc:778
lsst::afw::geom::Point2D GeomPoint
Definition: Wcs.cc:61
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
int status
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.

782  {
783  return skyToPixelImpl(sky1, sky2);
784 }
virtual geom::Point2D skyToPixelImpl(geom::Angle sky1, geom::Angle sky2) const
Definition: Wcs.cc:732
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.

771  {
773  return skyToPixelImpl(sky->getLongitude(), sky->getLatitude());
774 }
boost::shared_ptr< Coord > Ptr
Definition: Coord.h:72
afw::coord::Coord::Ptr convertCoordToSky(coord::Coord const &coord) const
Definition: Wcs.cc:778
virtual geom::Point2D skyToPixelImpl(geom::Angle sky1, geom::Angle sky2) const
Definition: Wcs.cc:732
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.

734  {
735  assert(_wcsInfo);
736 
737  double skyTmp[2];
738  double imgcrd[2];
739  double phi, theta;
740  double pixTmp[2];
741  /*
742  printf("_skyCoordsReversed: %c\n", (_skyCoordsReversed ? 'T' : 'F'));
743  printf("wcsinfo.lat: %i, lng: %i\n", _wcsInfo->lat, _wcsInfo->lng);
744  */
745  // WCSLib is smart enough to notice and handle crazy SDSS CTYPE1 = DEC--TAN,
746  // by recording the indices of the long and lat coordinates.
747  skyTmp[_wcsInfo->lng] = sky1.asDegrees();
748  skyTmp[_wcsInfo->lat] = sky2.asDegrees();
749 
750  int stat[1];
751  int status = 0;
752  status = wcss2p(_wcsInfo, 1, 2, skyTmp, &phi, &theta, imgcrd, pixTmp, stat);
753  if (status == 9) {
754  throw LSST_EXCEPT(except::DomainError,
755  (boost::format("sky coordinates %s, %s degrees is not valid for this WCS")
756  % sky1.asDegrees() % sky2.asDegrees()
757  ).str()
758  );
759  }
760  if (status > 0) {
761  throw LSST_EXCEPT(except::RuntimeError,
762  (boost::format("Error: wcslib returned a status code of %d at sky %s, %s deg: %s") %
763  status % sky1.asDegrees() % sky2.asDegrees() % wcs_errmsg[status]).str());
764  }
765 
766  // wcslib assumes 1-indexed coords
769 }
double asDegrees() const
Definition: Angle.h:124
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
const int fitsToLsstPixels
Definition: Wcs.cc:80
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
int status
const double PixelZeroPos
FITS uses 1.0, SDSS uses 0.5, LSST uses 0.0 (http://dev.lsstcorp.org/trac/wiki/BottomLeftPixelProposa...
Definition: ImageUtils.h:42
Point< double, 2 > Point2D
Definition: Point.h:286
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.

1078  {
1079  WcsPersistenceHelper const & keys = WcsPersistenceHelper::get();
1080  afw::table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
1081  PTR(afw::table::BaseRecord) record = catalog.addNew();
1082  record->set(keys.crval, getSkyOrigin()->getPosition(afw::geom::degrees));
1083  record->set(keys.crpix, getPixelOrigin());
1084  Eigen::Matrix2d cdIn = getCDMatrix();
1085  Eigen::Map<Eigen::Matrix2d> cdOut((*record)[keys.cd].getData());
1086  cdOut = cdIn;
1087  record->set(keys.ctype1, std::string(_wcsInfo[0].ctype[0]));
1088  record->set(keys.ctype2, std::string(_wcsInfo[0].ctype[1]));
1089  record->set(keys.equinox, _wcsInfo[0].equinox);
1090  record->set(keys.radesys, std::string(_wcsInfo[0].radesys));
1091  record->set(keys.cunit1, std::string(_wcsInfo[0].cunit[0]));
1092  record->set(keys.cunit2, std::string(_wcsInfo[0].cunit[1]));
1093  handle.saveCatalog(catalog);
1094 }
#define PTR(...)
Definition: base.h:41
table::Key< table::Array< double > > cd
Definition: Wcs.cc:1037
CatalogT< BaseRecord > BaseCatalog
Definition: fwd.h:61
table::PointKey< double > crval
Definition: Wcs.cc:1035
table::Key< std::string > radesys
Definition: Wcs.cc:1041
table::Key< std::string > ctype2
Definition: Wcs.cc:1039
AngleUnit const degrees
Definition: Angle.h:92
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
table::Key< double > equinox
Definition: Wcs.cc:1040
Eigen::Matrix2d getCDMatrix() const
Returns the CD matrix.
Definition: Wcs.cc:573
lsst::afw::geom::Point2D getPixelOrigin() const
Returns CRPIX (corrected to LSST convention).
Definition: Wcs.cc:565
lsst::afw::coord::Coord::Ptr getSkyOrigin() const
Returns CRVAL. This need not be the centre of the image.
Definition: Wcs.cc:560
table::Key< std::string > cunit1
Definition: Wcs.cc:1042
table::Key< std::string > ctype1
Definition: Wcs.cc:1038
table::Key< std::string > cunit2
Definition: Wcs.cc:1043
table::PointKey< double > crpix
Definition: Wcs.cc:1036

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: