31 #include "boost/format.hpp" 33 #include "wcslib/wcs.h" 34 #include "wcslib/wcsfix.h" 35 #include "wcslib/wcshdr.h" 68 void Wcs::_setWcslibParams() {
85 : daf::
base::Citizen(typeid(this)),
92 _coordSystem(coord::UNKNOWN)
119 if (ctype1[0] ==
'G') {
122 }
else if (ctype1[0] ==
'E') {
136 (
boost::format(
"Failed to setup wcs structure with wcsset. Status %d: %s") %
137 status % wcs_errmsg[status])
156 initWcsLib(crval, crpix, CD, ctype1, ctype2, equinox, raDecSys, cunits1, cunits2);
174 _hackHeader = _constHeader->deepCopy();
175 _constHeader = _hackHeader;
182 : _constHeader(header), _hackHeader() {}
189 HeaderAccess access(header);
197 if (access.toRead()->exists(key) && access.toRead()->typeOf(key) ==
typeid(
std::string)) {
198 auto equinoxStr = access.toRead()->getAsString(key).c_str();
199 double equinox = ::atof(equinoxStr[0] ==
'J' ? equinoxStr + 1 : equinoxStr);
200 access.toWrite()->set(key, equinox);
207 string msg =
"Could not parse FITS WCS: no header cards found";
214 if (!access.toRead()->exists(
"CRPIX1") && !access.toRead()->exists(
"CRPIX1a")) {
215 string msg =
"Neither CRPIX1 not CRPIX1a found";
219 if (!access.toRead()->exists(
"CRPIX2") && !access.toRead()->exists(
"CRPIX2a")) {
220 string msg =
"Neither CRPIX2 not CRPIX2a found";
225 if (!access.toRead()->exists(
"CRVAL1") && !access.toRead()->exists(
"CRVAL1a")) {
226 string msg =
"Neither CRVAL1 not CRVAL1a found";
230 if (!access.toRead()->exists(
"CRVAL2") && !access.toRead()->exists(
"CRVAL2a")) {
231 string msg =
"Neither CRVAL2 not CRVAL2a found";
242 if (access.toRead()->exists(
"CD1_1") || access.toRead()->exists(
"CD1_2") ||
243 access.toRead()->exists(
"CD2_1") || access.toRead()->exists(
"CD2_2")) {
244 for (
int i = 1; i <= 2; ++i) {
245 for (
int j = 1; j <= 2; ++j) {
247 if (access.toRead()->exists(key)) {
248 double const val = access.toRead()->getAsDouble(key);
249 access.toWrite()->remove(key);
250 access.toWrite()->add(
"X_" + key, val);
254 if (access.toRead()->exists(key)) {
255 double const val = access.toRead()->getAsDouble(key);
256 access.toWrite()->remove(key);
257 access.toWrite()->add(
"X_" + key, val);
267 char* hdrString =
const_cast<char*
>(metadataStr.
c_str());
273 if (pihStatus != 0) {
275 (
boost::format(
"Could not parse FITS WCS: wcspih status = %d (%s)") % pihStatus %
276 wcs_errmsg[pihStatus])
281 const int* naxes = NULL;
284 if (fixStatus != 0) {
286 errStream <<
"Could not parse FITS WCS: wcsfix failed " <<
std::endl;
287 for (
int ii = 0; ii < NWCSFIX; ++ii) {
288 if (stats[ii] >= 0) {
289 errStream <<
"\t" << ii <<
": " << stats[ii] <<
" " << wcsfix_errmsg[stats[ii]] <<
std::endl;
291 errStream <<
"\t" << ii <<
": " << stats[ii] <<
std::endl;
299 if (!(access.toRead()->exists(
"RADESYS") || access.toRead()->exists(
"RADESYSa"))) {
302 if (access.toRead()->exists(
"RADECSYS")) {
304 }
else if (access.toRead()->exists(
"EQUINOX") || access.toRead()->exists(
"EQUINOXa")) {
305 std::string const EQUINOX = access.toRead()->exists(
"EQUINOX") ?
"EQUINOX" :
"EQUINOXa";
306 double const equinox = access.toRead()->getAsDouble(EQUINOX);
307 if (equinox < 1984) {
330 double const* cdelt =
_wcsInfo->cdelt;
334 cd[0] = cdelt[0] * pc[0];
335 cd[1] = cdelt[0] * pc[1];
336 cd[2] = cdelt[1] * pc[2];
337 cd[3] = cdelt[1] * pc[3];
345 if ((CD.rows() != 2) || (CD.cols() != 2)) {
346 string msg =
"CD is not a 2x2 matrix";
351 bool isValid = (cunits1 ==
"deg");
352 isValid |= (cunits1 ==
"arcmin");
353 isValid |= (cunits1 ==
"arcsec");
354 isValid |= (cunits1 ==
"mas");
357 string msg =
"CUNITS1 must be one of {deg|arcmin|arcsec|mas}";
361 isValid = (cunits2 ==
"deg");
362 isValid |= (cunits2 ==
"arcmin");
363 isValid |= (cunits2 ==
"arcsec");
364 isValid |= (cunits2 ==
"mas");
367 string msg =
"CUNITS2 must be one of {deg|arcmin|arcsec|mas}";
372 _wcsInfo =
static_cast<struct wcsprm*
>(
malloc(
sizeof(
struct wcsprm)));
378 int status = wcsini(
true, 2,
_wcsInfo);
381 (
boost::format(
"Failed to allocate memory with wcsini. Status %d: %s") % status %
394 for (
int i = 0; i < 2; ++i) {
395 for (
int j = 0; j < 2; ++j) {
396 _wcsInfo->cd[(2 * i) + j] = CD(i, j);
425 (
boost::format(
"Failed to setup wcs structure with wcsset. Status %d: %s") %
426 status % wcs_errmsg[status])
448 for (
int i = 0; i != rhs.
_nWcsInfo; ++i) {
454 (
boost::format(
"Could not copy WCS: wcscopy status = %d : %s") % status %
469 if (&other ==
this)
return true;
480 inline bool compareArrays(
double const*
a,
double const*
b,
int n) {
481 for (
int i = 0; i < n; ++i)
482 if (a[i] != b[i])
return false;
486 template <
typename T>
487 inline bool compareStringArrays(T a, T b,
int n) {
488 for (
int i = 0; i < n; ++i)
489 if (
strcmp(a[i], b[i]) != 0)
return false;
493 #define CHECK_NULLS(a, b) \ 496 if ((b) == NULL) return true; \ 499 if ((b) == NULL) return false; \ 559 for (
int i = 0; i < naxis; ++i) {
560 for (
int j = 0; j < naxis; ++j) {
561 C(i, j) =
_wcsInfo->cd[(i * naxis) + j];
594 while (nQuarter < 0) {
611 switch (nQuarter % 4) {
619 _wcsInfo->crpix[0] = -crpy + dimensions.getY() + 1;
627 _wcsInfo->crpix[0] = -crpx + dimensions.getX() + 1;
628 _wcsInfo->crpix[1] = -crpy + dimensions.getY() + 1;
636 _wcsInfo->crpix[1] = -crpx + dimensions.getX() + 1;
664 static double square(
double x) {
return x *
x; }
672 const double side = 1;
691 double area =
sqrt(square(dx[1] * dy[2] - dx[2] * dy[1]) + square(dx[2] * dy[0] - dx[0] * dy[2]) +
692 square(dx[0] * dy[1] - dx[1] * dy[0]));
694 return area / square(side) * square(180. /
geom::PI);
717 status = wcss2p(
_wcsInfo, 1, 2, skyTmp, &phi, &theta, imgcrd, pixTmp, stat);
720 (
boost::format(
"sky coordinates %s, %s degrees is not valid for this WCS") %
726 (
boost::format(
"Error: wcslib returned a status code of %d at sky %s, %s deg: %s") %
762 skyTmp[
_wcsInfo->lng] = sky->getLongitude().asDegrees();
763 skyTmp[
_wcsInfo->lat] = sky->getLatitude().asDegrees();
768 imgcrd[0] = imgcrd[1] = -1e6;
773 status = wcss2p(
_wcsInfo, 1, 2, skyTmp, &phi, &theta, imgcrd, pixTmp, stat);
776 (
boost::format(
"Error: wcslib returned a status code of %d at sky %s, %s deg: %s") %
777 status % skyTmp[0] % skyTmp[1] % wcs_errmsg[status])
813 status = wcsp2s(
_wcsInfo, 1, 2, pixTmp, imgcrd, &phi, &theta, sky, &status);
816 (
boost::format(
"Error: wcslib returned a status code of %d at pixel %s, %s: %s") %
817 status % pixel1 % pixel2 % wcs_errmsg[status])
864 (
boost::format(
"Can't create Coord object: Unrecognised coordinate system %s") % coordSystem)
882 const double side = 10;
889 m(0, 0) = dsky10.first.asAngularUnits(skyUnit) / side;
890 m(0, 1) = dsky01.first.asAngularUnits(skyUnit) / side;
891 m(1, 0) = dsky10.second.asAngularUnits(skyUnit) / side;
892 m(1, 1) = dsky01.second.asAngularUnits(skyUnit) / side;
894 Eigen::Vector2d sky00v;
895 sky00v << sky00.getX(), sky00.getY();
896 Eigen::Vector2d pix00v;
897 pix00v << pix00.getX(), pix00.getY();
932 struct WcsPersistenceHelper {
936 table::Key<table::Array<double> >
cd;
944 static WcsPersistenceHelper
const&
get() {
945 static WcsPersistenceHelper instance;
950 WcsPersistenceHelper(
const WcsPersistenceHelper&) =
delete;
951 WcsPersistenceHelper&
operator=(
const WcsPersistenceHelper&) =
delete;
954 WcsPersistenceHelper(WcsPersistenceHelper&&) =
delete;
955 WcsPersistenceHelper&
operator=(WcsPersistenceHelper&&) =
delete;
958 WcsPersistenceHelper()
962 cd(schema.addField<table::Array<double> >(
963 "cd",
"linear transform matrix, ordered (1_1, 2_1, 1_2, 2_2)", 4)),
966 equinox(schema.addField<
double>(
"equinox",
"equinox of coordinates")),
967 radesys(schema.addField<
std::string>(
"radesys",
"coordinate system for equinox", 72)),
970 schema.getCitizen().markPersistent();
974 std::string getWcsPersistenceName() {
return "Wcs"; }
976 WcsFactory registration(getWcsPersistenceName());
985 WcsPersistenceHelper
const&
keys = WcsPersistenceHelper::get();
988 record->set(keys.crval.getX(),
_wcsInfo[0].crval[0]);
989 record->set(keys.crval.getY(),
_wcsInfo[0].crval[1]);
992 Eigen::Map<Eigen::Matrix2d> cdOut((*record)[keys.cd].getData());
996 record->set(keys.equinox,
_wcsInfo[0].equinox);
1004 if (
_wcsInfo[0].naxis != 2)
return false;
1029 WcsPersistenceHelper
const&
keys = WcsPersistenceHelper::get();
1034 Eigen::Matrix2d
cd = Eigen::Map<Eigen::Matrix2d const>(record[keys.cd].getData());
1036 record.
get(keys.ctype2), record.
get(keys.equinox), record.
get(keys.radesys),
1037 record.
get(keys.cunit1), record.
get(keys.cunit2));
1043 WcsPersistenceHelper
const&
keys = WcsPersistenceHelper::get();
1081 wcsMetaData->set(
"CRVAL1" + wcsName, x0,
"Column pixel of Reference Pixel");
1082 wcsMetaData->set(
"CRVAL2" + wcsName, y0,
"Row pixel of Reference Pixel");
1083 wcsMetaData->set(
"CRPIX1" + wcsName, 1,
"Column Pixel Coordinate of Reference");
1084 wcsMetaData->set(
"CRPIX2" + wcsName, 1,
"Row Pixel Coordinate of Reference");
1085 wcsMetaData->set(
"CTYPE1" + wcsName,
"LINEAR",
"Type of projection");
1086 wcsMetaData->set(
"CTYPE2" + wcsName,
"LINEAR",
"Type of projection");
1087 wcsMetaData->set(
"CUNIT1" + wcsName,
"PIXEL",
"Column unit");
1088 wcsMetaData->set(
"CUNIT2" + wcsName,
"PIXEL",
"Row unit");
1109 if (metadata->
getAsDouble(
"CRPIX1" + wcsName) == 1 &&
1111 x0 = metadata->
getAsInt(
"CRVAL1" + wcsName);
1112 y0 = metadata->
getAsInt(
"CRVAL2" + wcsName);
1116 metadata->
remove(
"CRVAL1" + wcsName);
1117 metadata->
remove(
"CRVAL2" + wcsName);
1118 metadata->
remove(
"CRPIX1" + wcsName);
1119 metadata->
remove(
"CRPIX2" + wcsName);
1120 metadata->
remove(
"CTYPE1" + wcsName);
1121 metadata->
remove(
"CTYPE1" + wcsName);
1122 metadata->
remove(
"CUNIT1" + wcsName);
1123 metadata->
remove(
"CUNIT2" + wcsName);
1156 metadata->remove(*ptr);
1168 : XYTransform(), _dst(dst), _src(src), _isSameSkySystem(dst->
isSameSkySystem(*src)) {}
1171 return std::make_shared<XYTransformFromWcsPair>(
_dst->clone(),
_src->clone());
1178 _src->pixelToSky(pixel[0], pixel[1], sky0, sky1);
1179 return _dst->skyToPixel(sky0, sky1);
1182 return _dst->skyToPixel(*coordPtr);
1190 _dst->pixelToSky(pixel[0], pixel[1], sky0, sky1);
1191 return _src->skyToPixel(sky0, sky1);
1194 return _src->skyToPixel(*coordPtr);
1200 return std::make_shared<XYTransformFromWcsPair>(
_src,
_dst);
double getAsDouble(std::string const &name) const
Get the last value for any arithmetic property name (possibly hierarchical).
std::shared_ptr< Coord > 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.
bool operator==(Wcs const &other) const
virtual bool isPersistable() const
Whether the Wcs is persistable using afw::table::io archives.
geom::Angle pixelScale() const
Returns the pixel scale [Angle/pixel].
std::shared_ptr< afw::coord::Coord > convertCoordToSky(coord::Coord const &coord) const
Given a Coord (as a shared pointer), return the sky position in the correct coordinate system for thi...
static PointKey addFields(Schema &schema, std::string const &name, std::string const &doc, std::string const &unit)
Add a pair of _x, _y fields to a Schema, and return a PointKey that points to them.
table::Key< std::string > name
afw::coord::CoordSystem getCoordSystem() const
WcsFactory(std::string const &name)
virtual geom::AffineTransform linearizeSkyToPixelInternal(geom::Point2D const &pix, coord::Coord const &coord, geom::AngleUnit skyUnit) const
Implementation for the overloaded public linearizeSkyToPixel methods, requiring both a pixel coordina...
An object passed to Persistable::write to allow it to persist itself.
virtual void remove(std::string const &name)
Remove all values for a property name (possibly hierarchical).
geom::Point2D skyToPixel(geom::Angle sky1, geom::Angle sky2) const
Convert from sky coordinates (e.g.
A custom container class for records, based on std::vector.
Class for storing ordered metadata with comments.
AngleUnit constexpr radians
constant with units of radians
bool _mayBePersistable() const
Perform basic checks on whether *this might be persistable.
A base class for factory classes used to reconstruct objects from records.
int getAsInt(std::string const &name) const
Get the last value for a bool/char/short/int property name (possibly hierarchical).
int stripWcsKeywords(std::shared_ptr< lsst::daf::base::PropertySet > const &metadata, std::shared_ptr< Wcs const > const &wcs)
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.
AngleUnit constexpr degrees
constant with units of degrees
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.
virtual std::shared_ptr< lsst::daf::base::PropertyList > getFitsMetadata() const
Return a PropertyList containing FITS header keywords that can be used to save the Wcs...
lsst::afw::geom::Point2D GeomPoint
table::PointKey< double > crval
bool isSameSkySystem(Wcs const &wcs) const
Return true if a WCS has the same coordinate system and equinox as this one.
geom::Point2I getImageXY0FromMetadata(daf::base::PropertySet &metadata, std::string const &wcsName, bool strip=false)
table::Key< std::string > radesys
std::pair< lsst::afw::geom::Angle, lsst::afw::geom::Angle > getTangentPlaneOffset(Coord const &c) const
Get the offset on the tangent plane.
Implementation of the WCS standard for a any projection.
lsst::daf::base::PropertyList PropertyList
table::Key< std::string > ctype2
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
A class used to convert scalar POD types such as double to Angle.
Rotation angle is unknown.
double pixArea(lsst::afw::geom::Point2D pix00) const
Sky area covered by a pixel at position pix00 in units of square degrees.
lsst::afw::geom::Point2D getPosition(lsst::afw::geom::AngleUnit unit=lsst::afw::geom::degrees) const
Return our contents in a Point2D object.
std::map< Citizen const *, CitizenInfo > table
double constexpr PI
The ratio of a circle's circumference to diameter.
lsst::daf::base::PropertySet PropertySet
bool isFlipped() const
Does the Wcs follow the convention of North=Up, East=Left?
afw::table::PointKey< int > dimensions
void initWcsLibFromFits(std::shared_ptr< lsst::daf::base::PropertySet const > const &fitsMetadata)
Parse a fits header, extract the relevant metadata and create a Wcs object.
A class representing an angle.
Wcs & operator=(const Wcs &)
A base class for image defects.
std::shared_ptr< coord::Coord > pixelToSky(double pix1, double pix2) const
Convert from pixel position to sky coordinates (e.g.
virtual bool _isSubset(Wcs const &other) const
geom::LinearTransform getLinearTransform() const
Return the linear part of the Wcs, the CD matrix in FITS-speak, as an AffineTransform.
Point< double, 2 > Point2D
std::shared_ptr< lsst::daf::base::PropertyList > createTrivialWcsAsPropertySet(std::string const &wcsName, int const x0=0, int const y0=0)
int contains(Schema const &other, int flags=EQUAL_KEYS) const
Test whether the given schema is a subset of this.
void shift(Extent< T, N > const &offset)
Shift the point by the given offset.
table::Key< double > equinox
const int fitsToLsstPixels
#define CHECK_NULLS(a, b)
virtual void shiftReferencePixel(double dx, double dy)
Move the pixel reference position by (dx, dy)
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Wcs()
Construct an invalid Wcs given no arguments.
Include files required for standard LSST Exception handling.
Schema getSchema() const
Return the Schema that holds this record's fields and keys.
lsst::afw::geom::Point2D getPixelOrigin() const
Returns CRPIX (corrected to LSST convention).
std::shared_ptr< afw::coord::Coord > 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...
void _initWcs()
Set some internal variables that we need to refer to.
geom::Point2D skyToIntermediateWorldCoord(coord::Coord const &coord) const
Convert from sky coordinates (e.g.
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
std::shared_ptr< lsst::afw::coord::Coord > getSkyOrigin() const
Returns CRVAL. This need not be the centre of the image.
virtual void rotateImageBy90(int nQuarter, lsst::afw::geom::Extent2I dimensions) const
Rotate image by nQuarter times 90 degrees.
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
virtual std::string getPythonModule() const
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
std::shared_ptr< Coord > convert(CoordSystem system, double epoch=2000) const
Convert to a specified Coord type at a specified epoch.
bool _skyAxesSwapped
if true then the sky axes are swapped
#define LSST_EXCEPT(type,...)
Create an exception with a given type and message and optionally other arguments (dependent on the ty...
A vector of catalogs used by Persistable.
Base class for all records.
int _wcshdrCtrl
Controls messages to stderr from wcshdr (0 for none); see wcshdr.h for details.
int _relax
Degree of permissiveness for wcspih (0 for strict); see wcshdr.h for details.
Citizen(const std::type_info &)
virtual void pixelToSkyImpl(double pixel1, double pixel2, geom::Angle skyTmp[2]) const
Worker routine for pixelToSky.
std::shared_ptr< lsst::afw::coord::Coord > CoordPtr
Class for storing generic metadata.
virtual geom::Point2D skyToPixelImpl(geom::Angle sky1, geom::Angle sky2) const
Worker routine for skyToPixel.
int _wcsfixCtrl
Do potentially unsafe translations of non-standard unit strings? 0/1 = no/yes.
virtual void write(OutputArchiveHandle &handle) const
Write the object to one or more catalogs.
table::Key< std::string > cunit1
CoordSystem makeCoordEnum(std::string const system)
A utility function to get the enum value of a coordinate system from a string name.
table::PointKey< int > pixel
virtual geom::AffineTransform linearizePixelToSkyInternal(geom::Point2D const &pix, coord::Coord const &coord, geom::AngleUnit skyUnit) const
Implementation for the overloaded public linearizePixelToSky methods, requiring both a pixel coordina...
table::Key< table::Array< double > > cd
virtual std::shared_ptr< table::io::Persistable > read(InputArchive const &archive, CatalogVector const &catalogs) const
Construct a new object from the given InputArchive and vector of catalogs.
double getEquinox() const
const int lsstToFitsPixels
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
Eigen::Matrix2d getCDMatrix() const
Returns the CD matrix.
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.
constexpr double asDegrees() const noexcept
Return an Angle's value in degrees.
This is the base class for spherical coordinates.
const double PixelZeroPos
position of center of pixel 0
virtual std::shared_ptr< Wcs > clone(void) const
table::Key< std::string > ctype1
coord::CoordSystem _coordSystem
virtual void flipImage(int flipLR, int flipTB, lsst::afw::geom::Extent2I dimensions) const
Flip CD matrix around the y-axis.
table::Key< std::string > cunit2
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
virtual std::string getPersistenceName() const
Return the unique name used to persist this object and look up its factory.
table::PointKey< double > crpix