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
Wcs.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
25 #include <iostream>
26 #include <sstream>
27 #include <cmath>
28 #include <cstring>
29 
30 #include "boost/format.hpp"
31 
32 #include "wcslib/wcs.h"
33 #include "wcslib/wcsfix.h"
34 #include "wcslib/wcshdr.h"
35 
36 #include "lsst/daf/base.h"
37 #include "lsst/daf/base/Citizen.h"
40 #include "lsst/pex/exceptions.h"
42 #include "lsst/afw/image/Wcs.h"
43 #include "lsst/afw/coord/Coord.h"
44 #include "lsst/afw/geom/Angle.h"
49 
50 namespace except = lsst::pex::exceptions;
51 namespace afwImg = lsst::afw::image;
52 namespace afwCoord = lsst::afw::coord;
53 namespace afwGeom = lsst::afw::geom;
54 
55 
56 using namespace std;
57 
64 
65 //The amount of space allocated to strings in wcslib
66 const int STRLEN = 72;
67 
68 //Set internal params for wcslib
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 }
78 
79 const int lsstToFitsPixels = +1;
80 const int fitsToLsstPixels = -1;
81 
82 //
83 // Constructors
84 //
85 
86 
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 }
95 
96 
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 }
114 
115 /*
116  * Set some internal variables that we need to refer to
117  */
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 }
134 
149 Wcs::Wcs(GeomPoint const & crval, GeomPoint const & crpix, Eigen::Matrix2d const & CD,
150  std::string const & ctype1, std::string const & ctype2,
151  double equinox, std::string const & raDecSys,
152  std::string const & cunits1, std::string const & cunits2
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 }
170 
171 
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 }
350 
361 void Wcs::initWcsLib(GeomPoint const & crval, GeomPoint const & crpix, Eigen::Matrix2d const & CD,
362  std::string const & ctype1, std::string const & ctype2,
363  double equinox, std::string const & raDecSys,
364  std::string const & cunits1, std::string const & cunits2) {
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 }
451 
452 
454 Wcs::Wcs(afwImg::Wcs const & rhs) :
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 }
485 
486 bool Wcs::operator==(Wcs const & other) const {
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 }
494 
495 // convenience functions and a macro for implementing isSubset
496 namespace {
497 
498 inline bool compareArrays(double const * a, double const * b, int n) {
499  for (int i = 0; i < n; ++i) if (a[i] != b[i]) return false;
500  return true;
501 }
502 
503 template <typename T>
504 inline bool compareStringArrays(T a, T b, int n) {
505  for (int i = 0; i < n; ++i) if (strcmp(a[i], b[i]) != 0) return false;
506  return true;
507 }
508 
509 #define CHECK_NULLS(a, b) \
510  do { \
511  if ((a) == NULL) { \
512  if ((b) == NULL) return true; \
513  return false; \
514  } \
515  if ((b) == NULL) return false; \
516  } while (false)
517 
518 } // anonymous
519 
520 bool Wcs::_isSubset(Wcs const & rhs) const {
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 }
544 
546  if (_wcsInfo != NULL) {
547  wcsvfree(&_nWcsInfo, &_wcsInfo);
548  }
549 }
550 
551 
552 Wcs::Ptr Wcs::clone(void) const {
553  return Wcs::Ptr(new Wcs(*this));
554 }
555 
556 //
557 // Accessors
558 //
559 
561  assert(_wcsInfo);
562  return makeCorrectCoord(_wcsInfo->crval[0] * afwGeom::degrees, _wcsInfo->crval[1] * afwGeom::degrees);
563 }
564 
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 }
572 
573 Eigen::Matrix2d Wcs::getCDMatrix() const {
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 }
590 
591 void Wcs::flipImage(int flipLR, int flipTB, afwGeom::Extent2I dimensions) const {
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 }
613 
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 }
667 
670 }
671 
672 bool Wcs::isFlipped() const {
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 }
687 
688 static double square(double x) {
689  return x*x;
690 }
691 
692 double Wcs::pixArea(GeomPoint pix00
693  ) const {
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 }
724 
726  return sqrt(pixArea(getPixelOrigin())) * afwGeom::degrees;
727 }
728 
729 /*
730  * Worker routine for skyToPixel
731  */
732 GeomPoint Wcs::skyToPixelImpl(afwGeom::Angle sky1, // RA (or, more generally, longitude)
733  afwGeom::Angle sky2 // Dec (or latitude)
734  ) const {
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 }
770 
773  return skyToPixelImpl(sky->getLongitude(), sky->getLatitude());
774 }
775 
776 
779  return coord.convert(_coordSystem);
780 }
781 
783  return skyToPixelImpl(sky1, sky2);
784 }
785 
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 }
826 
827 /*
828  * Worker routine for pixelToSky
829  */
830 void
831 Wcs::pixelToSkyImpl(double pixel1, double pixel2, afwGeom::Angle skyTmp[2]) const
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 }
853 
855  return pixelToSky(pixel.getX(), pixel.getY());
856 }
857 
858 CoordPtr Wcs::pixelToSky(double pixel1, double pixel2) const {
859  assert(_wcsInfo);
860 
861  afwGeom::Angle skyTmp[2];
862  pixelToSkyImpl(pixel1, pixel2, skyTmp);
863  return makeCorrectCoord(skyTmp[0], skyTmp[1]);
864 }
865 
866 void Wcs::pixelToSky(double pixel1, double pixel2, afwGeom::Angle& sky1, afwGeom::Angle& sky2) const {
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 }
876 
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 }
932 
933 
935  lsst::afw::coord::Coord const & coord,
937 ) const {
938  return linearizePixelToSkyInternal(skyToPixel(coord), coord, skyUnit);
939 }
941  GeomPoint const & pix,
943 ) const {
944  return linearizePixelToSkyInternal(pix, *pixelToSky(pix), skyUnit);
945 }
946 
947 /*
948  * Implementation for the overloaded public linearizePixelToSky methods, requiring both a pixel coordinate
949  * and the corresponding sky coordinate.
950  */
952  GeomPoint const & pix00,
953  lsst::afw::coord::Coord const & coord,
955 ) const {
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 }
980 
982  lsst::afw::coord::Coord const & coord,
984 ) const {
985  return linearizeSkyToPixelInternal(skyToPixel(coord), coord, skyUnit);
986 }
987 
989  GeomPoint const & pix,
991 ) const {
992  return linearizeSkyToPixelInternal(pix, *pixelToSky(pix), skyUnit);
993 }
994 
1000  GeomPoint const & pix00,
1001  lsst::afw::coord::Coord const & coord,
1003 ) const {
1004  lsst::afw::geom::AffineTransform inverse = linearizePixelToSkyInternal(pix00, coord, skyUnit);
1005  return inverse.invert();
1006 }
1007 
1009 {
1011 }
1012 
1013 // -------- table-based persistence -------------------------------------------------------------------------
1014 
1015 namespace lsst { namespace afw { namespace image {
1016 
1018 public:
1019 
1020  explicit WcsFactory(std::string const & name) : table::io::PersistableFactory(name) {}
1021 
1022  virtual PTR(table::io::Persistable) read(
1023  InputArchive const & archive,
1024  CatalogVector const & catalogs
1025  ) const;
1026 
1027 };
1028 
1029 namespace {
1030 
1031 // Read-only singleton struct containing the schema and keys that a simple Wcs is mapped
1032 // to in record persistence.
1033 struct WcsPersistenceHelper : private boost::noncopyable {
1034  table::Schema schema;
1035  table::PointKey<double> crval;
1036  table::PointKey<double> crpix;
1037  table::Key< table::Array<double> > cd;
1038  table::Key<std::string> ctype1;
1039  table::Key<std::string> ctype2;
1040  table::Key<double> equinox;
1041  table::Key<std::string> radesys;
1042  table::Key<std::string> cunit1;
1043  table::Key<std::string> cunit2;
1044 
1045  static WcsPersistenceHelper const & get() {
1046  static WcsPersistenceHelper instance;
1047  return instance;
1048  };
1049 
1050 private:
1051  WcsPersistenceHelper() :
1052  schema(),
1053  crval(table::PointKey<double>::addFields(schema, "crval", "celestial reference point", "degrees")),
1054  crpix(table::PointKey<double>::addFields(schema, "crpix", "pixel reference point", "pixels")),
1055  cd(schema.addField< table::Array<double> >(
1056  "cd", "linear transform matrix, ordered (1_1, 2_1, 1_2, 2_2)", 4)),
1057  ctype1(schema.addField< std::string >("ctype1", "coordinate type", 72)),
1058  ctype2(schema.addField< std::string >("ctype2", "coordinate type", 72)),
1059  equinox(schema.addField< double >("equinox", "equinox of coordinates")),
1060  radesys(schema.addField< std::string >("radesys", "coordinate system for equinox", 72)),
1061  cunit1(schema.addField< std::string >("cunit1", "coordinate units", 72)),
1062  cunit2(schema.addField< std::string >("cunit2", "coordinate units", 72))
1063  {
1064  schema.getCitizen().markPersistent();
1065  }
1066 };
1067 
1068 std::string getWcsPersistenceName() { return "Wcs"; }
1069 
1070 WcsFactory registration(getWcsPersistenceName());
1071 
1072 } // anonymous
1073 
1074 std::string Wcs::getPersistenceName() const { return getWcsPersistenceName(); }
1075 
1076 std::string Wcs::getPythonModule() const { return "lsst.afw.image"; }
1077 
1078 void Wcs::write(OutputArchiveHandle & handle) const {
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 }
1095 
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 }
1103 
1104 bool Wcs::isPersistable() const {
1105  if (!_mayBePersistable()) {
1106  return false;
1107  }
1108  // The current table persistence only works for TAN and TAN-SIP projections
1109  return false;
1110 }
1111 
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 }
1139 
1141 WcsFactory::read(InputArchive const & inputs, CatalogVector const & catalogs) const {
1142  WcsPersistenceHelper const & keys = WcsPersistenceHelper::get();
1143  LSST_ARCHIVE_ASSERT(catalogs.size() >= 1u);
1144  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
1145  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
1146  PTR(Wcs) result(new Wcs(catalogs.front().front()));
1147  return result;
1148 }
1149 
1150 }}} // namespace lsst::afw::image
1151 
1152 // ----------------------------------------------------------------------------------------------------------
1153 
1154 //
1155 //Mutators
1156 //
1157 
1158 
1164 void Wcs::shiftReferencePixel(double dx, double dy) {
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 }
1172 
1173 /************************************************************************************************************/
1174 /*
1175  * Now WCSA, pixel coordinates, but allowing for X0 and Y0
1176  */
1177 namespace lsst {
1178 namespace afw {
1179 namespace image {
1180 
1181 namespace detail {
1182 
1187 createTrivialWcsAsPropertySet(std::string const& wcsName,
1188  int const x0,
1189  int const y0
1190  ) {
1192 
1193  wcsMetaData->set("CRVAL1" + wcsName, x0, "Column pixel of Reference Pixel");
1194  wcsMetaData->set("CRVAL2" + wcsName, y0, "Row pixel of Reference Pixel");
1195  wcsMetaData->set("CRPIX1" + wcsName, 1, "Column Pixel Coordinate of Reference");
1196  wcsMetaData->set("CRPIX2" + wcsName, 1, "Row Pixel Coordinate of Reference");
1197  wcsMetaData->set("CTYPE1" + wcsName, "LINEAR", "Type of projection");
1198  wcsMetaData->set("CTYPE2" + wcsName, "LINEAR", "Type of projection");
1199  wcsMetaData->set("CUNIT1" + wcsName, "PIXEL", "Column unit");
1200  wcsMetaData->set("CUNIT2" + wcsName, "PIXEL", "Row unit");
1201 
1202  return wcsMetaData;
1203 }
1210 afwGeom::Point2I getImageXY0FromMetadata(std::string const& wcsName,
1211  lsst::daf::base::PropertySet *metadata
1212  ) {
1213 
1214  int x0 = 0; // Our value of X0
1215  int y0 = 0; // Our value of Y0
1216 
1217  try {
1218  //
1219  // Only use WCS if CRPIX[12] == 1 and CRVAL[12] is present
1220  //
1221  if (metadata->getAsDouble("CRPIX1" + wcsName) == 1 &&
1222  metadata->getAsDouble("CRPIX2" + wcsName) == 1) {
1223 
1224  x0 = metadata->getAsInt("CRVAL1" + wcsName);
1225  y0 = metadata->getAsInt("CRVAL2" + wcsName);
1226  //
1227  // OK, we've got it. Remove it from the header
1228  //
1229  metadata->remove("CRVAL1" + wcsName);
1230  metadata->remove("CRVAL2" + wcsName);
1231  metadata->remove("CRPIX1" + wcsName);
1232  metadata->remove("CRPIX2" + wcsName);
1233  metadata->remove("CTYPE1" + wcsName);
1234  metadata->remove("CTYPE1" + wcsName);
1235  metadata->remove("CUNIT1" + wcsName);
1236  metadata->remove("CUNIT2" + wcsName);
1237  }
1238  } catch(lsst::pex::exceptions::NotFoundError &) {
1239  ; // OK, not present
1240  }
1241 
1242  return afwGeom::Point2I(x0, y0);
1243 }
1244 
1254  CONST_PTR(Wcs) const& wcs
1255  )
1256 {
1257  PTR(lsst::daf::base::PropertySet) wcsMetadata = wcs->getFitsMetadata();
1258  std::vector<std::string> paramNames = wcsMetadata->paramNames();
1259  paramNames.push_back("CDELT1");
1260  paramNames.push_back("CDELT2");
1261  paramNames.push_back("LTV1");
1262  paramNames.push_back("LTV2");
1263  paramNames.push_back("PC001001");
1264  paramNames.push_back("PC001002");
1265  paramNames.push_back("PC002001");
1266  paramNames.push_back("PC002002");
1267  for (std::vector<std::string>::const_iterator ptr = paramNames.begin(); ptr != paramNames.end(); ++ptr) {
1268  metadata->remove(*ptr);
1269  }
1270 
1271  return 0; // would be ncard if remove returned a status
1272 }
1273 
1274 
1275 }}}}
1276 
1277 
1278 
1279 // -------------------------------------------------------------------------------------------------
1280 //
1281 // XYTransformFromWcsPair
1282 
1283 
1285  : XYTransform(), _dst(dst), _src(src)
1286 { }
1287 
1288 
1290 {
1291  return boost::make_shared<XYTransformFromWcsPair>(_dst->clone(), _src->clone());
1292 }
1293 
1294 
1296 {
1297  //
1298  // TODO there is an alternate version of pixelToSky() which is designated for the
1299  // "knowledgeable user in need of performance". This is probably better, but first I need
1300  // to understand exactly which checks are needed (e.g. I think we need to check by hand
1301  // that both Wcs's use the same celestial coordinate system)
1302  //
1303  PTR(afw::coord::Coord) x = _src->pixelToSky(pixel);
1304  return _dst->skyToPixel(*x);
1305 }
1306 
1308 {
1309  PTR(afw::coord::Coord) x = _dst->pixelToSky(pixel);
1310  return _src->skyToPixel(*x);
1311 }
1312 
1314 {
1315  // just swap src, dst
1316  return boost::make_shared<XYTransformFromWcsPair> (_src, _dst);
1317 }
int _wcsfixCtrl
Do potentially unsafe translations of non-standard unit strings? 0/1 = no/yes.
Definition: Wcs.h:394
geom::Point2I getImageXY0FromMetadata(std::string const &wcsName, lsst::daf::base::PropertySet *metadata)
Definition: Wcs.cc:1210
const int lsstToFitsPixels
Definition: Wcs.cc:79
Interface for WcsFormatter class.
boost::shared_ptr< Coord > Ptr
Definition: Coord.h:72
Schema getSchema() const
Return the Schema that holds this record&#39;s fields and keys.
Definition: BaseRecord.h:48
virtual std::string getPersistenceName() const
Return the unique name used to persist this object and look up its factory.
Definition: Wcs.cc:1074
boost::shared_ptr< Wcs const > _src
Definition: Wcs.h:447
#define PTR(...)
Definition: base.h:41
table::Key< std::string > name
Definition: ApCorrMap.cc:71
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
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.
Definition: Wcs.cc:934
WcsFactory(std::string const &name)
Definition: Wcs.cc:1020
table::Key< table::Array< double > > cd
Definition: Wcs.cc:1037
An object passed to Persistable::write to allow it to persist itself.
virtual void remove(std::string const &name)
Definition: PropertySet.cc:754
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:47
AngleUnit const radians
constant with units of radians
Definition: Angle.h:91
AffineTransform const invert() const
virtual void pixelToSkyImpl(double pixel1, double pixel2, geom::Angle skyTmp[2]) const
Definition: Wcs.cc:831
coord::CoordSystem _coordSystem
Definition: Wcs.h:397
void _setWcslibParams()
Definition: Wcs.cc:69
A custom container class for records, based on std::vector.
Definition: Catalog.h:94
Class for storing ordered metadata with comments.
Definition: PropertyList.h:81
boost::shared_ptr< Wcs > Ptr
Definition: Wcs.h:113
A base class for factory classes used to reconstruct objects from records.
Definition: Persistable.h:231
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
#define CONST_PTR(...)
Definition: base.h:47
afw::coord::Coord::Ptr convertCoordToSky(coord::Coord const &coord) const
Definition: Wcs.cc:778
std::pair< lsst::afw::geom::Angle, lsst::afw::geom::Angle > getTangentPlaneOffset(Coord const &c) const
Get the offset on the tangent plane.
Definition: Coord.cc:748
virtual bool _isSubset(Wcs const &other) const
Definition: Wcs.cc:520
XYTransformFromWcsPair: An XYTransform obtained by putting two Wcs objects &quot;back to back&quot;...
Definition: Wcs.h:432
lsst::afw::geom::Point2D GeomPoint
Definition: Wcs.cc:61
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
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
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
tbl::Key< int > wcs
static lsst::daf::base::PropertyList::Ptr generatePropertySet(lsst::afw::image::Wcs const &wcs)
virtual geom::AffineTransform linearizeSkyToPixelInternal(geom::Point2D const &pix, coord::Coord const &coord, geom::AngleUnit skyUnit) const
Definition: Wcs.cc:999
table::Key< std::string > radesys
Definition: Wcs.cc:1041
geom::Point2D skyToIntermediateWorldCoord(coord::Coord const &coord) const
Convert from sky coordinates (e.g. RA/dec) to intermediate world coordinates.
Definition: Wcs.cc:786
virtual ~Wcs()
Definition: Wcs.cc:545
virtual Ptr clone(void) const
Definition: Wcs.cc:552
Implementation of the WCS standard for a any projection.
Definition: Wcs.h:107
Point< int, 2 > Point2I
Definition: PSF.h:39
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::daf::base::PropertyList PropertyList
Definition: Wcs.cc:59
table::Key< std::string > ctype2
Definition: Wcs.cc:1039
A class used to convert scalar POD types such as double to Angle.
Definition: Angle.h:71
int const x0
Definition: saturated.cc:45
Image utility functions.
double asDegrees() const
Definition: Angle.h:124
virtual void shiftReferencePixel(double dx, double dy)
Move the pixel reference position by (dx, dy)
Definition: Wcs.cc:1164
boost::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Definition: Catalog.h:469
virtual geom::Point2D skyToPixelImpl(geom::Angle sky1, geom::Angle sky2) const
Definition: Wcs.cc:732
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
A base class for objects that can be persisted via afw::table::io Archive classes.
Definition: Persistable.h:74
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition: BaseRecord.h:123
A 2D linear coordinate transformation.
lsst::daf::base::PropertySet PropertySet
Definition: Wcs.cc:58
virtual std::string getPythonModule() const
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
Definition: Wcs.cc:1076
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...
Definition: Wcs.cc:668
virtual Point2D reverseTransform(Point2D const &pixel) const
Definition: Wcs.cc:1307
int getAsInt(std::string const &name) const
Definition: PropertySet.cc:333
afw::table::PointKey< int > dimensions
AngleUnit const degrees
Definition: Angle.h:92
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
int contains(Schema const &other, int flags=EQUAL_KEYS) const
Test whether the given schema is a subset of this.
virtual bool isPersistable() const
Whether the Wcs is persistable using afw::table::io archives.
Definition: Wcs.cc:1104
boost::shared_ptr< PropertyList > Ptr
Definition: PropertyList.h:84
An affine coordinate transformation consisting of a linear transformation and an offset.
virtual void rotateImageBy90(int nQuarter, lsst::afw::geom::Extent2I dimensions) const
Rotate image by nQuarter times 90 degrees.
Definition: Wcs.cc:614
table::Key< double > equinox
Definition: Wcs.cc:1040
bool operator==(Wcs const &other) const
Definition: Wcs.cc:486
#define CHECK_NULLS(a, b)
Definition: Wcs.cc:509
Eigen::Matrix2d getCDMatrix() const
Returns the CD matrix.
Definition: Wcs.cc:573
geom::Angle pixelScale() const
Returns the pixel scale [Angle/pixel].
Definition: Wcs.cc:725
const int STRLEN
Definition: Wcs.cc:66
lsst::afw::geom::Point2D getPixelOrigin() const
Returns CRPIX (corrected to LSST convention).
Definition: Wcs.cc:565
const int fitsToLsstPixels
Definition: Wcs.cc:80
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
virtual void flipImage(int flipLR, int flipTB, lsst::afw::geom::Extent2I dimensions) const
Flip CD matrix around the y-axis.
Definition: Wcs.cc:591
boost::shared_ptr< lsst::daf::base::PropertyList > createTrivialWcsAsPropertySet(std::string const &wcsName, int const x0=0, int const y0=0)
Definition: Wcs.cc:1187
lsst::afw::coord::Coord::Ptr getSkyOrigin() const
Returns CRVAL. This need not be the centre of the image.
Definition: Wcs.cc:560
virtual geom::AffineTransform linearizePixelToSkyInternal(geom::Point2D const &pix, coord::Coord const &coord, geom::AngleUnit skyUnit) const
Definition: Wcs.cc:951
tbl::Schema schema
void shift(Extent< T, N > const &offset)
Shift the point by the given offset.
Definition: Point.h:113
int x
lsst::afw::image::Wcs Wcs
Definition: Wcs.cc:60
void _initWcs()
Definition: Wcs.cc:118
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
A vector of catalogs used by Persistable.
Definition: CatalogVector.h:26
int _relax
Degree of permissiveness for wcspih (0 for strict); see wcshdr.h for details.
Definition: Wcs.h:393
bool isFlipped() const
Definition: Wcs.cc:672
Base class for all records.
Definition: BaseRecord.h:27
tuple m
Definition: lsstimport.py:48
geom::LinearTransform getLinearTransform() const
Definition: Wcs.cc:1008
Wcs()
Construct an invalid Wcs given no arguments.
Definition: Wcs.cc:88
int status
double getAsDouble(std::string const &name) const
Definition: PropertySet.cc:406
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
Class for storing generic metadata.
Definition: PropertySet.h:82
Virtual base class for 2D transforms.
Definition: XYTransform.h:48
XYTransformFromWcsPair(boost::shared_ptr< Wcs const > dst, boost::shared_ptr< Wcs const > src)
Definition: Wcs.cc:1284
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
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.
Definition: Wcs.cc:981
afw::table::Key< double > b
lsst::afw::image::XYTransformFromWcsPair XYTransformFromWcsPair
Definition: Wcs.cc:63
table::Key< std::string > cunit1
Definition: Wcs.cc:1042
Point< double, 2 > Point2D
Definition: Point.h:286
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
table::PointKey< int > pixel
std::string formatFitsProperties(boost::shared_ptr< lsst::daf::base::PropertySet const > const &prop)
double const PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:18
virtual void write(OutputArchiveHandle &handle) const
Write the object to one or more catalogs.
Definition: Wcs.cc:1078
A multi-catalog archive object used to load table::io::Persistable objects.
Definition: InputArchive.h:28
Coord::Ptr convert(CoordSystem system) const
Convert to a specified Coord type.
Definition: Coord.cc:637
int stripWcsKeywords(boost::shared_ptr< lsst::daf::base::PropertySet > const &metadata, boost::shared_ptr< Wcs const > const &wcs)
Definition: Wcs.cc:1253
int const y0
Definition: saturated.cc:45
int _wcshdrCtrl
Controls messages to stderr from wcshdr (0 for none); see wcshdr.h for details.
Definition: Wcs.h:395
lsst::afw::coord::Coord::Ptr CoordPtr
Definition: Wcs.cc:62
ImageT val
Definition: CR.cc:154
virtual Point2D forwardTransform(Point2D const &pixel) const
virtuals for forward and reverse transforms
Definition: Wcs.cc:1295
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
boost::shared_ptr< Wcs const > _dst
Definition: Wcs.h:446
int countFitsHeaderCards(boost::shared_ptr< lsst::daf::base::PropertySet const > const &prop)
table::Key< std::string > ctype1
Definition: Wcs.cc:1038
lsst::afw::geom::Point2D getPosition(lsst::afw::geom::AngleUnit unit=lsst::afw::geom::degrees) const
Return our contents in a Point2D object.
Definition: Coord.cc:477
bool _mayBePersistable() const
Perform basic checks on whether *this might be persistable.
Definition: Wcs.cc:1096
Include files required for standard LSST Exception handling.
table::Key< std::string > cunit2
Definition: Wcs.cc:1043
table::PointKey< double > crpix
Definition: Wcs.cc:1036
Functions to handle coordinates.
Formatting utilities.