LSSTApplications  11.0-22-g33de520,13.0+153,14.0+52,14.0+57,14.0-1-g013352c+36,14.0-1-g13ef843+9,14.0-1-g4b114ac+14,14.0-1-g7257b6a+12,14.0-1-g8b7e855+51,14.0-13-g7a60b79+2,14.0-14-g87d16e8+10,14.0-14-gbf7a6f8a,14.0-17-g4f4ea82+5,14.0-2-g319577b+11,14.0-2-ga5af9b6+10,14.0-22-gc48c03f+3,14.0-3-g20413be+3,14.0-46-g76222d5f+3,14.0-47-g0a51fac97,14.0-5-g744ff5f+2,14.0-5-g86eb1bd+31,14.0-6-gd5b81a9+6,14.0-6-ge2c9487+42,14.0-8-g7f6dd6b+6,14.0-8-gb81b6e9+4,14.0-9-g11010eb,14.0-9-g330837b+5
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 #include <limits>
30 
31 #include "boost/format.hpp"
32 
33 #include "wcslib/wcs.h"
34 #include "wcslib/wcsfix.h"
35 #include "wcslib/wcshdr.h"
36 
37 #include "lsst/daf/base.h"
38 #include "lsst/daf/base/Citizen.h"
41 #include "lsst/pex/exceptions.h"
43 #include "lsst/afw/image/Wcs.h"
44 #include "lsst/afw/coord/Coord.h"
45 #include "lsst/afw/geom/Angle.h"
50 
51 namespace except = lsst::pex::exceptions;
52 
53 using namespace std;
54 
59 
60 // The amount of space allocated to strings in wcslib
61 const int STRLEN = 72;
62 
63 namespace lsst {
64 namespace afw {
65 namespace image {
66 
67 // Set internal params for wcslib
68 void Wcs::_setWcslibParams() {
69  _wcsfixCtrl = // ctrl for wcsfix
70  2; // Translate "H" to "h"
71  _wcshdrCtrl = // ctrl for wcspih
72  2; // Report each rejected keyrecord and the reason why it was rejected
73  _relax = // relax parameter for wcspih;
74  WCSHDR_all; // Accept all extensions recognized by the parser
75 }
76 
77 const int lsstToFitsPixels = +1;
78 const int fitsToLsstPixels = -1;
79 
80 //
81 // Constructors
82 //
83 
84 Wcs::Wcs()
85  : daf::base::Citizen(typeid(this)),
86  _wcsInfo(NULL),
87  _nWcsInfo(0),
88  _relax(0),
89  _wcsfixCtrl(0),
90  _wcshdrCtrl(0),
91  _nReject(0),
92  _coordSystem(coord::UNKNOWN) // set by _initWcs
93 {
95  _initWcs();
96 }
97 
99  : daf::base::Citizen(typeid(this)),
100  _wcsInfo(NULL),
101  _nWcsInfo(0),
102  _relax(0),
103  _wcsfixCtrl(0),
104  _wcshdrCtrl(0),
105  _nReject(0),
106  _coordSystem(coord::UNKNOWN) // set by _initWcs
107 {
109 
110  initWcsLibFromFits(fitsMetadata);
111  _initWcs();
112 }
113 
115  // first four characters of CTYPE1 (name of first axis)
116  std::string ctype1 = std::string(_wcsInfo->ctype[0]).substr(0, 4);
117 
118  if (_wcsInfo) {
119  if (ctype1[0] == 'G') {
121  _skyAxesSwapped = (ctype1[2] == 'A'); // GLAT instead of GLON
122  } else if (ctype1[0] == 'E') {
124  _skyAxesSwapped = (ctype1[2] == 'A'); // ELAT instead of ELON
125  } else {
127  _skyAxesSwapped = (ctype1[0] == 'D'); // DEC instead of RA
128  }
129 
130  // tell WCSlib that values have been updated
131  _wcsInfo->flag = 0;
132  // and then tell it to do its internal magic.
133  int status = wcsset(_wcsInfo);
134  if (status != 0) {
136  (boost::format("Failed to setup wcs structure with wcsset. Status %d: %s") %
137  status % wcs_errmsg[status])
138  .str());
139  }
140  }
141 }
142 
143 Wcs::Wcs(GeomPoint const& crval, GeomPoint const& crpix, Eigen::Matrix2d const& CD, std::string const& ctype1,
144  std::string const& ctype2, double equinox, std::string const& raDecSys, std::string const& cunits1,
145  std::string const& cunits2)
146  : daf::base::Citizen(typeid(this)),
147  _wcsInfo(NULL),
148  _nWcsInfo(0),
149  _relax(0),
150  _wcsfixCtrl(0),
151  _wcshdrCtrl(0),
152  _nReject(0),
153  _coordSystem(coord::UNKNOWN) // set by _initWcs
154 {
156  initWcsLib(crval, crpix, CD, ctype1, ctype2, equinox, raDecSys, cunits1, cunits2);
157  _initWcs();
158 }
159 
161  /*
162  * Access control for the input header
163  *
164  * We want to hack up the input, and in order to do so we need to do a deep copy on it.
165  * We only want to do that copy once, and would like to avoid doing it altogether.
166  */
167  class HeaderAccess {
168  public:
169  // Return a readable version of the metadata
170  std::shared_ptr<daf::base::PropertySet const> const& toRead() { return _constHeader; }
171  // Return a writable version of the metadata
172  std::shared_ptr<daf::base::PropertySet> const& toWrite() {
173  if (!_hackHeader) {
174  _hackHeader = _constHeader->deepCopy();
175  _constHeader = _hackHeader;
176  }
177  return _hackHeader;
178  }
179 
180  // Ctor
181  HeaderAccess(std::shared_ptr<daf::base::PropertySet const> const& header)
182  : _constHeader(header), _hackHeader() {}
183 
184  private:
187  };
188 
189  HeaderAccess access(header);
190 
191  // Some headers (e.g. SDSS ones from FNAL) have EQUINOX as a string. Fix this,
192  // as wcslib 4.4.4 refuses to handle it. Furthermore, some headers (e.g. Skymapper)
193  // use a string but prepend 'J': "J2000.0" -- if we don't handle this, we deduce
194  // an equinox of 0.0
195  {
196  std::string const& key = "EQUINOX";
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);
201  }
202  }
203 
204  //Check header isn't empty
205  int nCards = lsst::afw::formatters::countFitsHeaderCards(*(access.toRead()));
206  if (nCards <= 0) {
207  string msg = "Could not parse FITS WCS: no header cards found";
209  }
210 
211  // While the standard does not insist on CRVAL and CRPIX being present, it
212  // is almost certain their absence indicates a problem.
213  // Check for CRPIX
214  if (!access.toRead()->exists("CRPIX1") && !access.toRead()->exists("CRPIX1a")) {
215  string msg = "Neither CRPIX1 not CRPIX1a found";
217  }
218 
219  if (!access.toRead()->exists("CRPIX2") && !access.toRead()->exists("CRPIX2a")) {
220  string msg = "Neither CRPIX2 not CRPIX2a found";
222  }
223 
224  // And the same for CRVAL
225  if (!access.toRead()->exists("CRVAL1") && !access.toRead()->exists("CRVAL1a")) {
226  string msg = "Neither CRVAL1 not CRVAL1a found";
228  }
229 
230  if (!access.toRead()->exists("CRVAL2") && !access.toRead()->exists("CRVAL2a")) {
231  string msg = "Neither CRVAL2 not CRVAL2a found";
233  }
234  /*
235  * According to Greisen and Calabretta (A&A 395, 1061–1075 (2002)) it's illegal to mix PCi_j and CDi_j
236  * headers; unfortunately Subaru puts both in its headers. It actually uses PC001002 instead of PC1_2
237  * (dating to a proposed FITS standard from 1996) and at least sometimes fails to include CDELT[12],
238  * so the CD and PC matrices are inconsistent
239  *
240  * If we detect any part of a CD matrix, delete all PC matrices
241  */
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) {
246  std::string key = (boost::format("PC%i_%i") % j % i).str();
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);
251  }
252 
253  key = (boost::format("PC%03d%03d") % j % i).str();
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);
258  }
259  }
260  }
261  }
262 
263  //Pass the header into wcslib's formatter to extract & setup the Wcs. First need
264  //to convert to a C style string, so the compile doesn't complain about constness
265  std::string metadataStr = formatters::formatFitsProperties(*(access.toRead()));
266  // We own the data, and wcslib is slack about constness, so no qualms with casting away const
267  char* hdrString = const_cast<char*>(metadataStr.c_str());
268  // printf("wcspih string:\n%s\n", hdrString);
269 
270  nCards = formatters::countFitsHeaderCards(*(access.toRead())); // we may have dropped some
271  int pihStatus = wcspih(hdrString, nCards, _relax, _wcshdrCtrl, &_nReject, &_nWcsInfo, &_wcsInfo);
272 
273  if (pihStatus != 0) {
275  (boost::format("Could not parse FITS WCS: wcspih status = %d (%s)") % pihStatus %
276  wcs_errmsg[pihStatus])
277  .str());
278  }
279 
280  // Run wcsfix on _wcsInfo to try and fix any problems it knows about.
281  const int* naxes = NULL; // should be {NAXIS1, NAXIS2, ...} to check cylindrical projections
282  int stats[NWCSFIX]; // status returns from wcsfix
283  int fixStatus = wcsfix(_wcsfixCtrl, naxes, _wcsInfo, stats);
284  if (fixStatus != 0) {
285  std::stringstream errStream;
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;
290  } else {
291  errStream << "\t" << ii << ": " << stats[ii] << std::endl;
292  }
293  }
294  }
295 
296  // The Wcs standard requires a default value for RADESYS if the keyword
297  // doesn't exist in header, but wcslib doesn't set it. So we do so here. This code
298  // conforms to Calabretta & Greisen 2002 \S 3.1
299  if (!(access.toRead()->exists("RADESYS") || access.toRead()->exists("RADESYSa"))) {
300  // If RADECSYS exists, use that (counter to Calabretta & Greisen 2002 \S 3.1, but commonly used).
301  // If equinox exist and < 1984, use FK4. If >= 1984, use FK5
302  if (access.toRead()->exists("RADECSYS")) {
303  strncpy(_wcsInfo->radesys, access.toRead()->getAsString("RADECSYS").c_str(), STRLEN);
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) {
308  strncpy(_wcsInfo->radesys, "FK4", STRLEN);
309  } else {
310  strncpy(_wcsInfo->radesys, "FK5", STRLEN);
311  }
312  } else {
313  // If Equinox doesn't exist, default to ICRS
314  strncpy(_wcsInfo->radesys, "ICRS", STRLEN);
315  }
316  }
317  // strip trailing whitespace
318  {
319  for (int i = strlen(_wcsInfo->radesys) - 1; i >= 0; i--) {
320  if (isspace(_wcsInfo->radesys[i])) {
321  _wcsInfo->radesys[i] = '\0';
322  }
323  }
324  }
325  //
326  // If there are no CDi_j cards in the header, set CDi_j from PCi_j
327  // CDi_j == CDELTi*PCi_j
328  //
329  if ((_wcsInfo->altlin & 2) == 0) { // no CDi_j cards were found in the header
330  double const* cdelt = _wcsInfo->cdelt;
331  double const* pc = _wcsInfo->pc;
332  double* cd = _wcsInfo->cd;
333 
334  cd[0] = cdelt[0] * pc[0]; // 1_1
335  cd[1] = cdelt[0] * pc[1]; // 1_2
336  cd[2] = cdelt[1] * pc[2]; // 2_1
337  cd[3] = cdelt[1] * pc[3]; // 2_2
338  }
339 }
340 
341 void Wcs::initWcsLib(GeomPoint const& crval, GeomPoint const& crpix, Eigen::Matrix2d const& CD,
342  std::string const& ctype1, std::string const& ctype2, double equinox,
343  std::string const& raDecSys, std::string const& cunits1, std::string const& cunits2) {
344  // Check CD is a valid size
345  if ((CD.rows() != 2) || (CD.cols() != 2)) {
346  string msg = "CD is not a 2x2 matrix";
348  }
349 
350  // Check that cunits are legitimate values
351  bool isValid = (cunits1 == "deg");
352  isValid |= (cunits1 == "arcmin");
353  isValid |= (cunits1 == "arcsec");
354  isValid |= (cunits1 == "mas");
355 
356  if (!isValid) {
357  string msg = "CUNITS1 must be one of {deg|arcmin|arcsec|mas}";
359  }
360 
361  isValid = (cunits2 == "deg");
362  isValid |= (cunits2 == "arcmin");
363  isValid |= (cunits2 == "arcsec");
364  isValid |= (cunits2 == "mas");
365 
366  if (!isValid) {
367  string msg = "CUNITS2 must be one of {deg|arcmin|arcsec|mas}";
369  }
370 
371  // Initialise the wcs struct
372  _wcsInfo = static_cast<struct wcsprm*>(malloc(sizeof(struct wcsprm)));
373  if (_wcsInfo == NULL) {
374  throw LSST_EXCEPT(except::MemoryError, "Cannot allocate WCS info");
375  }
376 
377  _wcsInfo->flag = -1;
378  int status = wcsini(true, 2, _wcsInfo); // 2 indicates a naxis==2, a two dimensional image
379  if (status != 0) {
381  (boost::format("Failed to allocate memory with wcsini. Status %d: %s") % status %
382  wcs_errmsg[status])
383  .str());
384  }
385 
386  // Set crval, crpix and CD. Internally to the class, we use fits units for consistency with
387  // wcslib.
388  _wcsInfo->crval[0] = crval.getX();
389  _wcsInfo->crval[1] = crval.getY();
390  _wcsInfo->crpix[0] = crpix.getX() + lsstToFitsPixels;
391  _wcsInfo->crpix[1] = crpix.getY() + lsstToFitsPixels;
392 
393  // Set the CD matrix
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);
397  }
398  }
399 
400  // Specify that we have a CD matrix, but no PC or CROTA
401  _wcsInfo->altlin = 2;
402  _wcsInfo->flag = 0; // values have been updated
403 
404  // This is a work around for what I think is a bug in wcslib. ->types is neither
405  // initialised or set to NULL by default, so if I try to delete a Wcs object,
406  // wcslib then attempts to free non-existent space, and the code can crash.
407  _wcsInfo->types = NULL;
408 
409  // Set the coordinate system
410  strncpy(_wcsInfo->ctype[0], ctype1.c_str(), STRLEN);
411  strncpy(_wcsInfo->ctype[1], ctype2.c_str(), STRLEN);
412  strncpy(_wcsInfo->radesys, raDecSys.c_str(), STRLEN);
413  _wcsInfo->equinox = equinox;
414 
415  // Set the units
416  strncpy(_wcsInfo->cunit[0], cunits1.c_str(), STRLEN);
417  strncpy(_wcsInfo->cunit[1], cunits2.c_str(), STRLEN);
418 
419  _nWcsInfo = 1; // Specify that we have only one coordinate representation
420 
421  // Tell wcslib that we are need to set up internal values
422  status = wcsset(_wcsInfo);
423  if (status != 0) {
425  (boost::format("Failed to setup wcs structure with wcsset. Status %d: %s") %
426  status % wcs_errmsg[status])
427  .str());
428  }
429 }
430 
431 Wcs::Wcs(Wcs const& rhs)
432  : daf::base::Citizen(typeid(this)),
433  _wcsInfo(NULL),
434  _nWcsInfo(rhs._nWcsInfo),
435  _relax(rhs._relax),
438  _nReject(rhs._nReject),
439  _coordSystem(coord::UNKNOWN) // set by _initWcs
440 {
441  if (rhs._nWcsInfo > 0) {
442  _wcsInfo = static_cast<struct wcsprm*>(calloc(rhs._nWcsInfo, sizeof(struct wcsprm)));
443  if (_wcsInfo == NULL) {
444  throw LSST_EXCEPT(pex::exceptions::MemoryError, "Cannot allocate WCS info");
445  }
446 
447  int alloc = 1; // Unconditionally allocate memory when calling
448  for (int i = 0; i != rhs._nWcsInfo; ++i) {
449  _wcsInfo[i].flag = -1;
450  int status = wcscopy(alloc, &rhs._wcsInfo[i], &_wcsInfo[i]);
451  if (status != 0) {
452  wcsvfree(&i, &_wcsInfo);
454  (boost::format("Could not copy WCS: wcscopy status = %d : %s") % status %
455  wcs_errmsg[status])
456  .str());
457  }
458  }
459  }
460  _initWcs();
461 }
462 // Delegate to copy-constructor for backwards compatibility
463 Wcs::Wcs(Wcs&& rhs) : Wcs(rhs) {}
464 
465 Wcs& Wcs::operator=(const Wcs&) = default;
466 Wcs& Wcs::operator=(Wcs&&) = default;
467 
468 bool Wcs::operator==(Wcs const& other) const {
469  if (&other == this) return true;
470  // We do a bidirectional test with a virtual member function in case one of us is a derived
471  // class with members we don't know about here.
472  // This is not the most efficient possible implementation, but I think it's the easiest one
473  // with which to ensure correctness, and I think that's more important in this case.
474  return this->_isSubset(other) && other._isSubset(*this);
475 }
476 
477 // convenience functions and a macro for implementing isSubset
478 namespace {
479 
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;
483  return true;
484 }
485 
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;
490  return true;
491 }
492 
493 #define CHECK_NULLS(a, b) \
494  do { \
495  if ((a) == NULL) { \
496  if ((b) == NULL) return true; \
497  return false; \
498  } \
499  if ((b) == NULL) return false; \
500  } while (false)
501 
502 } // anonymous
503 
504 bool Wcs::_isSubset(Wcs const& rhs) const {
506  CHECK_NULLS(_wcsInfo->crval, rhs._wcsInfo->crval);
507  CHECK_NULLS(_wcsInfo->cd, rhs._wcsInfo->cd);
508  CHECK_NULLS(_wcsInfo->crpix, rhs._wcsInfo->crpix);
509  CHECK_NULLS(_wcsInfo->cunit, rhs._wcsInfo->cunit);
510  CHECK_NULLS(_wcsInfo->ctype, rhs._wcsInfo->ctype);
511  return _nWcsInfo == rhs._nWcsInfo && _coordSystem == rhs._coordSystem &&
512  _wcsInfo->naxis == rhs._wcsInfo->naxis && _wcsInfo->equinox == rhs._wcsInfo->equinox &&
513  _wcsInfo->altlin == rhs._wcsInfo->altlin &&
514  compareArrays(_wcsInfo->crval, rhs._wcsInfo->crval, 2) &&
515  compareArrays(_wcsInfo->crpix, rhs._wcsInfo->crpix, 2) &&
516  compareArrays(_wcsInfo->cd, rhs._wcsInfo->cd, 4) &&
517  compareStringArrays(_wcsInfo->cunit, rhs._wcsInfo->cunit, 2) &&
518  compareStringArrays(_wcsInfo->ctype, rhs._wcsInfo->ctype, 2) &&
519  skyToPixel(_wcsInfo->crval[0] * geom::degrees, _wcsInfo->crval[1] * geom::degrees) ==
520  rhs.skyToPixel(_wcsInfo->crval[0] * geom::degrees, _wcsInfo->crval[1] * geom::degrees) &&
521  *pixelToSky(_wcsInfo->crpix[0], _wcsInfo->crpix[1]) ==
522  *rhs.pixelToSky(_wcsInfo->crpix[0], _wcsInfo->crpix[1]);
523 }
524 
526  if (_wcsInfo != NULL) {
527  wcsvfree(&_nWcsInfo, &_wcsInfo);
528  }
529 }
530 
531 std::shared_ptr<Wcs> Wcs::clone(void) const { return std::shared_ptr<Wcs>(new Wcs(*this)); }
532 
533 //
534 // Accessors
535 //
536 
538  assert(_wcsInfo);
539  return makeCorrectCoord(_wcsInfo->crval[0] * geom::degrees, _wcsInfo->crval[1] * geom::degrees);
540 }
541 
543  assert(_wcsInfo);
544  // Convert from fits units back to lsst units
545  double p1 = _wcsInfo->crpix[0] + fitsToLsstPixels;
546  double p2 = _wcsInfo->crpix[1] + fitsToLsstPixels;
547  return geom::Point2D(p1, p2);
548 }
549 
550 Eigen::Matrix2d Wcs::getCDMatrix() const {
551  assert(_wcsInfo);
552  int const naxis = _wcsInfo->naxis;
553 
554  // If naxis != 2, I'm not sure if any of what follows is correct
555  assert(naxis == 2);
556 
557  Eigen::Matrix2d C;
558 
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];
562  }
563  }
564 
565  return C;
566 }
567 
568 void Wcs::flipImage(int flipLR, int flipTB, geom::Extent2I dimensions) const {
569  assert(_wcsInfo);
570 
571  int const naxis = _wcsInfo->naxis;
572 
573  // If naxis != 2, I'm not sure if any of what follows is correct
574  assert(naxis == 2);
575  // Origin is at (1,1). Adjust to avoid off by one.
576  if (flipLR) {
577  _wcsInfo->cd[0] = -_wcsInfo->cd[0];
578  _wcsInfo->cd[2] = -_wcsInfo->cd[2];
579  _wcsInfo->crpix[0] = -_wcsInfo->crpix[0] + dimensions.getX() + 1;
580  }
581  if (flipTB) {
582  _wcsInfo->cd[1] = -_wcsInfo->cd[1];
583  _wcsInfo->cd[3] = -_wcsInfo->cd[3];
584  _wcsInfo->crpix[1] = -_wcsInfo->crpix[1] + dimensions.getY() + 1;
585  }
586 
587  // tells libwcs to invalidate cached data, since transformation has been modified
588  _wcsInfo->flag = 0;
589 }
590 
591 void Wcs::rotateImageBy90(int nQuarter, geom::Extent2I dimensions) const {
592  assert(_wcsInfo);
593 
594  while (nQuarter < 0) {
595  nQuarter += 4;
596  }
597 
598  int const naxis = _wcsInfo->naxis;
599 
600  // If naxis != 2, I'm not sure if any of what follows is correct
601  assert(naxis == 2);
602  double a = _wcsInfo->cd[0];
603  double b = _wcsInfo->cd[1];
604  double c = _wcsInfo->cd[2];
605  double d = _wcsInfo->cd[3];
606  double crpx = _wcsInfo->crpix[0];
607  double crpy = _wcsInfo->crpix[1];
608  // Origin is at (1,1). Adjust to avoid off by one.
609  // E.g. CRPIX one pixel off the UR of the grid should go to
610  //(0,0) for nQuarter=2
611  switch (nQuarter % 4) {
612  case 0:
613  break;
614  case 1:
615  _wcsInfo->cd[0] = -b;
616  _wcsInfo->cd[1] = a;
617  _wcsInfo->cd[2] = -d;
618  _wcsInfo->cd[3] = c;
619  _wcsInfo->crpix[0] = -crpy + dimensions.getY() + 1;
620  _wcsInfo->crpix[1] = crpx;
621  break;
622  case 2:
623  _wcsInfo->cd[0] = -a;
624  _wcsInfo->cd[1] = -b;
625  _wcsInfo->cd[2] = -c;
626  _wcsInfo->cd[3] = -d;
627  _wcsInfo->crpix[0] = -crpx + dimensions.getX() + 1;
628  _wcsInfo->crpix[1] = -crpy + dimensions.getY() + 1;
629  break;
630  case 3:
631  _wcsInfo->cd[0] = b;
632  _wcsInfo->cd[1] = -a;
633  _wcsInfo->cd[2] = d;
634  _wcsInfo->cd[3] = -c;
635  _wcsInfo->crpix[0] = crpy;
636  _wcsInfo->crpix[1] = -crpx + dimensions.getX() + 1;
637  break;
638  }
639 
640  // tells libwcs to invalidate cached data, since transformation has been modified
641  _wcsInfo->flag = 0;
642 }
643 
646 }
647 
648 bool Wcs::isFlipped() const {
649  // We calculate the determinant of the CD (i.e the rotation and scaling)
650  // matrix. If this determinant is positive, then the image can be rotated
651  // to a position where increasing the right ascension and declination
652  // increases the horizontal and vertical pixel position. In this case the
653  // image is flipped.
654  assert(_wcsInfo);
655  double det = (_wcsInfo->cd[0] * _wcsInfo->cd[3]) - (_wcsInfo->cd[1] * _wcsInfo->cd[2]);
656 
657  if (det == 0) {
658  throw(LSST_EXCEPT(except::RuntimeError, "Wcs CD matrix is singular"));
659  }
660 
661  return (det > 0);
662 }
663 
664 static double square(double x) { return x * x; }
665 
666 double Wcs::pixArea(GeomPoint pix00) const {
667  //
668  // Figure out the (0, 0), (0, 1), and (1, 0) ra/dec coordinates of the corners of a square drawn in pixel
669  // It'd be better to centre the square at sky00, but that would involve another conversion between sky and
670  // pixel coordinates so I didn't bother
671  //
672  const double side = 1; // length of the square's sides in pixels
673 
674  // Work in 3-space to avoid RA wrapping and pole issues.
675  geom::Point3D v0 = pixelToSky(pix00)->getVector();
676 
677  // Step by "side" in x and y pixel directions...
678  GeomPoint px(pix00);
679  GeomPoint py(pix00);
680  px.shift(geom::Extent2D(side, 0));
681  py.shift(geom::Extent2D(0, side));
682  // Push the points through the WCS, and find difference in 3-space.
683  geom::Extent3D dx = pixelToSky(px)->getVector() - v0;
684  geom::Extent3D dy = pixelToSky(py)->getVector() - v0;
685 
686  // Compute |cross product| = area of parallelogram with sides dx,dy
687  // FIXME -- this is slightly incorrect -- it's making the small-angle
688  // approximation, taking the distance *through* the unit sphere
689  // rather than over its surface.
690  // This is in units of ~radians^2
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]));
693 
694  return area / square(side) * square(180. / geom::PI);
695 }
696 
698 
700  assert(_wcsInfo);
701 
702  double skyTmp[2];
703  double imgcrd[2];
704  double phi, theta;
705  double pixTmp[2];
706  /*
707  printf("_skyCoordsReversed: %c\n", (_skyCoordsReversed ? 'T' : 'F'));
708  printf("wcsinfo.lat: %i, lng: %i\n", _wcsInfo->lat, _wcsInfo->lng);
709  */
710  // WCSLib is smart enough to notice and handle crazy SDSS CTYPE1 = DEC--TAN,
711  // by recording the indices of the long and lat coordinates.
712  skyTmp[_wcsInfo->lng] = sky1.asDegrees();
713  skyTmp[_wcsInfo->lat] = sky2.asDegrees();
714 
715  int stat[1];
716  int status = 0;
717  status = wcss2p(_wcsInfo, 1, 2, skyTmp, &phi, &theta, imgcrd, pixTmp, stat);
718  if (status == 9) {
720  (boost::format("sky coordinates %s, %s degrees is not valid for this WCS") %
721  sky1.asDegrees() % sky2.asDegrees())
722  .str());
723  }
724  if (status > 0) {
726  (boost::format("Error: wcslib returned a status code of %d at sky %s, %s deg: %s") %
727  status % sky1.asDegrees() % sky2.asDegrees() % wcs_errmsg[status])
728  .str());
729  }
730 
731  // wcslib assumes 1-indexed coords
732  return geom::Point2D(pixTmp[0] + PixelZeroPos + fitsToLsstPixels,
733  pixTmp[1] + PixelZeroPos + fitsToLsstPixels);
734 }
735 
738  return skyToPixelImpl(sky->getLongitude(), sky->getLatitude());
739 }
740 
742  return coord.convert(_coordSystem, _wcsInfo->equinox);
743 }
744 
745 GeomPoint Wcs::skyToPixel(geom::Angle sky1, geom::Angle sky2) const { return skyToPixelImpl(sky1, sky2); }
746 
748  assert(_wcsInfo);
749 
751  double skyTmp[2];
752  double imgcrd[2];
753  double phi, theta;
754  double pixTmp[2];
755 
756  /*
757  printf("skyToIWC: _coordSystem = %i\n", (int)_coordSystem);
758  printf("coord (%.3f, %.3f)\n", coord->getLongitude().asDegrees(), coord->getLatitude().asDegrees());
759  printf("->sky (%.3f, %.3f)\n", sky->getLongitude().asDegrees(), sky->getLatitude().asDegrees());
760  */
761 
762  skyTmp[_wcsInfo->lng] = sky->getLongitude().asDegrees();
763  skyTmp[_wcsInfo->lat] = sky->getLatitude().asDegrees();
764 
765  // Estimate pixel coordinates
766  int stat[1];
767  int status = 0;
768  imgcrd[0] = imgcrd[1] = -1e6;
769  /*
770  printf(" skyTmp[] = (%.3f, %.3f)\n", skyTmp[0], skyTmp[1]);
771  printf(" _wcsInfo->lng,lat = %i, %i\n", _wcsInfo->lng, _wcsInfo->lat);
772  */
773  status = wcss2p(_wcsInfo, 1, 2, skyTmp, &phi, &theta, imgcrd, pixTmp, stat);
774  if (status > 0) {
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])
778  .str());
779  }
780  /*
781  printf("->iwc (%.3f, %.3f)\n", imgcrd[0], imgcrd[1]);
782  printf("-> pix (%.2f, %.2f)\n", pixTmp[0], pixTmp[1]);
783  std::shared_ptr<coord::Coord> crval = getSkyOrigin();
784  printf("(crval is (%.3f, %.3f))\n", crval->getLongitude().asDegrees(), crval->getLatitude().asDegrees());
785  */
786  return GeomPoint(imgcrd[0], imgcrd[1]);
787 }
788 
789 double Wcs::getEquinox() const {
790  if (_wcsInfo == nullptr) {
792  }
793  return _wcsInfo->equinox;
794 }
795 
796 bool Wcs::isSameSkySystem(Wcs const& wcs) const {
797  if (_isIcrs() && wcs._isIcrs()) {
798  return true;
799  }
800  return (getCoordSystem() == wcs.getCoordSystem()) && (getEquinox() == wcs.getEquinox());
801 }
802 
803 void Wcs::pixelToSkyImpl(double pixel1, double pixel2, geom::Angle skyTmp[2]) const {
804  assert(_wcsInfo);
805 
806  // wcslib assumes 1-indexed coordinates
807  double pixTmp[2] = {pixel1 - PixelZeroPos + lsstToFitsPixels, pixel2 - PixelZeroPos + lsstToFitsPixels};
808  double imgcrd[2];
809  double phi, theta;
810 
811  double sky[2];
812  int status = 0;
813  status = wcsp2s(_wcsInfo, 1, 2, pixTmp, imgcrd, &phi, &theta, sky, &status);
814  if (status > 0) {
816  (boost::format("Error: wcslib returned a status code of %d at pixel %s, %s: %s") %
817  status % pixel1 % pixel2 % wcs_errmsg[status])
818  .str());
819  }
820  // FIXME -- _wcsInfo.lat, _wcsInfo.lng ?
821  skyTmp[0] = sky[0] * geom::degrees;
822  skyTmp[1] = sky[1] * geom::degrees;
823 }
824 
825 CoordPtr Wcs::pixelToSky(GeomPoint const& pixel) const { return pixelToSky(pixel.getX(), pixel.getY()); }
826 
827 CoordPtr Wcs::pixelToSky(double pixel1, double pixel2) const {
828  assert(_wcsInfo);
829 
830  geom::Angle skyTmp[2];
831  pixelToSkyImpl(pixel1, pixel2, skyTmp);
832  return makeCorrectCoord(skyTmp[0], skyTmp[1]);
833 }
834 
835 void Wcs::pixelToSky(double pixel1, double pixel2, geom::Angle& sky1, geom::Angle& sky2) const {
836  geom::Angle skyTmp[2];
837  // HACK -- we shouldn't need to initialize these -- pixelToSkyImpl() sets them unless an
838  // exception is thrown -- but be safe.
839  skyTmp[0] = 0. * geom::radians;
840  skyTmp[1] = 0. * geom::radians;
841  pixelToSkyImpl(pixel1, pixel2, skyTmp);
842  sky1 = skyTmp[0];
843  sky2 = skyTmp[1];
844 }
845 
847  auto coordSystem = getCoordSystem();
848  if ((coordSystem == coord::ICRS) || (coordSystem == coord::GALACTIC)) {
849  // equinox not relevant
850  if (_skyAxesSwapped) {
851  return coord::makeCoord(coordSystem, sky1, sky0);
852  } else {
853  return coord::makeCoord(coordSystem, sky0, sky1);
854  }
855  } else if ((coordSystem == coord::FK5) || (coordSystem == coord::ECLIPTIC)) {
856  if (_skyAxesSwapped) {
857  return coord::makeCoord(coordSystem, sky1, sky0, _wcsInfo->equinox);
858  } else {
859  return coord::makeCoord(coordSystem, sky0, sky1, _wcsInfo->equinox);
860  }
861  }
862  throw LSST_EXCEPT(
864  (boost::format("Can't create Coord object: Unrecognised coordinate system %s") % coordSystem)
865  .str());
866 }
867 
869  return linearizePixelToSkyInternal(skyToPixel(coord), coord, skyUnit);
870 }
872  return linearizePixelToSkyInternal(pix, *pixelToSky(pix), skyUnit);
873 }
874 
876  geom::AngleUnit skyUnit) const {
877  //
878  // Figure out the (0, 0), (0, 1), and (1, 0) ra/dec coordinates of the corners of a square drawn in pixel
879  // It'd be better to centre the square at sky00, but that would involve another conversion between sky and
880  // pixel coordinates so I didn't bother
881  //
882  const double side = 10; // length of the square's sides in pixels
883  GeomPoint const sky00 = coord.getPosition(skyUnit);
884  typedef std::pair<geom::Angle, geom::Angle> AngleAngle;
885  AngleAngle const dsky10 = coord.getTangentPlaneOffset(*pixelToSky(pix00 + geom::Extent2D(side, 0)));
886  AngleAngle const dsky01 = coord.getTangentPlaneOffset(*pixelToSky(pix00 + geom::Extent2D(0, side)));
887 
888  Eigen::Matrix2d m;
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;
893 
894  Eigen::Vector2d sky00v;
895  sky00v << sky00.getX(), sky00.getY();
896  Eigen::Vector2d pix00v;
897  pix00v << pix00.getX(), pix00.getY();
898  // return geom::AffineTransform(m, geom::Extent2D(sky00v - m * pix00v));
899  return geom::AffineTransform(m, (sky00v - m * pix00v));
900 }
901 
903  return linearizeSkyToPixelInternal(skyToPixel(coord), coord, skyUnit);
904 }
905 
907  return linearizeSkyToPixelInternal(pix, *pixelToSky(pix), skyUnit);
908 }
909 
911  geom::AngleUnit skyUnit) const {
912  geom::AffineTransform inverse = linearizePixelToSkyInternal(pix00, coord, skyUnit);
913  return inverse.invert();
914 }
915 
917 
918 // -------- table-based persistence -------------------------------------------------------------------------
919 
921 public:
922  explicit WcsFactory(std::string const& name) : table::io::PersistableFactory(name) {}
923 
924  virtual std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
925  CatalogVector const& catalogs) const;
926 };
927 
928 namespace {
929 
930 // Read-only singleton struct containing the schema and keys that a simple Wcs is mapped
931 // to in record persistence.
932 struct WcsPersistenceHelper {
933  table::Schema schema;
934  table::PointKey<double> crval;
935  table::PointKey<double> crpix;
936  table::Key<table::Array<double> > cd;
937  table::Key<std::string> ctype1;
938  table::Key<std::string> ctype2;
939  table::Key<double> equinox;
940  table::Key<std::string> radesys;
941  table::Key<std::string> cunit1;
942  table::Key<std::string> cunit2;
943 
944  static WcsPersistenceHelper const& get() {
945  static WcsPersistenceHelper instance;
946  return instance;
947  };
948 
949  // No copying
950  WcsPersistenceHelper(const WcsPersistenceHelper&) = delete;
951  WcsPersistenceHelper& operator=(const WcsPersistenceHelper&) = delete;
952 
953  // No moving
954  WcsPersistenceHelper(WcsPersistenceHelper&&) = delete;
955  WcsPersistenceHelper& operator=(WcsPersistenceHelper&&) = delete;
956 
957 private:
958  WcsPersistenceHelper()
959  : schema(),
960  crval(table::PointKey<double>::addFields(schema, "crval", "celestial reference point", "deg")),
961  crpix(table::PointKey<double>::addFields(schema, "crpix", "pixel reference point", "pixel")),
962  cd(schema.addField<table::Array<double> >(
963  "cd", "linear transform matrix, ordered (1_1, 2_1, 1_2, 2_2)", 4)),
964  ctype1(schema.addField<std::string>("ctype1", "coordinate type", 72)),
965  ctype2(schema.addField<std::string>("ctype2", "coordinate type", 72)),
966  equinox(schema.addField<double>("equinox", "equinox of coordinates")),
967  radesys(schema.addField<std::string>("radesys", "coordinate system for equinox", 72)),
968  cunit1(schema.addField<std::string>("cunit1", "coordinate units", 72)),
969  cunit2(schema.addField<std::string>("cunit2", "coordinate units", 72)) {
970  schema.getCitizen().markPersistent();
971  }
972 };
973 
974 std::string getWcsPersistenceName() { return "Wcs"; }
975 
976 WcsFactory registration(getWcsPersistenceName());
977 
978 } // anonymous
979 
980 std::string Wcs::getPersistenceName() const { return getWcsPersistenceName(); }
981 
982 std::string Wcs::getPythonModule() const { return "lsst.afw.image"; }
983 
984 void Wcs::write(OutputArchiveHandle& handle) const {
985  WcsPersistenceHelper const& keys = WcsPersistenceHelper::get();
986  afw::table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
988  record->set(keys.crval.getX(), _wcsInfo[0].crval[0]);
989  record->set(keys.crval.getY(), _wcsInfo[0].crval[1]);
990  record->set(keys.crpix, getPixelOrigin());
991  Eigen::Matrix2d cdIn = getCDMatrix();
992  Eigen::Map<Eigen::Matrix2d> cdOut((*record)[keys.cd].getData());
993  cdOut = cdIn;
994  record->set(keys.ctype1, std::string(_wcsInfo[0].ctype[0]));
995  record->set(keys.ctype2, std::string(_wcsInfo[0].ctype[1]));
996  record->set(keys.equinox, _wcsInfo[0].equinox);
997  record->set(keys.radesys, std::string(_wcsInfo[0].radesys));
998  record->set(keys.cunit1, std::string(_wcsInfo[0].cunit[0]));
999  record->set(keys.cunit2, std::string(_wcsInfo[0].cunit[1]));
1000  handle.saveCatalog(catalog);
1001 }
1002 
1004  if (_wcsInfo[0].naxis != 2) return false;
1005  if (std::strcmp(_wcsInfo[0].cunit[0], "deg") != 0) return false;
1006  if (std::strcmp(_wcsInfo[0].cunit[1], "deg") != 0) return false;
1007 
1008  return true;
1009 }
1010 
1011 bool Wcs::isPersistable() const {
1012  if (!_mayBePersistable()) {
1013  return false;
1014  }
1015  // The current table persistence only works for TAN and TAN-SIP projections
1016  return false;
1017 }
1018 
1020  : daf::base::Citizen(typeid(this)),
1021  _wcsInfo(NULL),
1022  _nWcsInfo(0),
1023  _relax(0),
1024  _wcsfixCtrl(0),
1025  _wcshdrCtrl(0),
1026  _nReject(0),
1027  _coordSystem(static_cast<afw::coord::CoordSystem>(-1)) // set by _initWcs
1028 {
1029  WcsPersistenceHelper const& keys = WcsPersistenceHelper::get();
1030  if (!record.getSchema().contains(keys.schema)) {
1031  throw LSST_EXCEPT(afw::table::io::MalformedArchiveError, "Incorrect schema for Wcs persistence");
1032  }
1033  _setWcslibParams();
1034  Eigen::Matrix2d cd = Eigen::Map<Eigen::Matrix2d const>(record[keys.cd].getData());
1035  initWcsLib(record.get(keys.crval), record.get(keys.crpix), cd, record.get(keys.ctype1),
1036  record.get(keys.ctype2), record.get(keys.equinox), record.get(keys.radesys),
1037  record.get(keys.cunit1), record.get(keys.cunit2));
1038  _initWcs();
1039 }
1040 
1042  CatalogVector const& catalogs) const {
1043  WcsPersistenceHelper const& keys = WcsPersistenceHelper::get();
1044  LSST_ARCHIVE_ASSERT(catalogs.size() >= 1u);
1045  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
1046  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
1047  std::shared_ptr<Wcs> result(new Wcs(catalogs.front().front()));
1048  return result;
1049 }
1050 
1051 // ----------------------------------------------------------------------------------------------------------
1052 
1053 //
1054 // Mutators
1055 //
1056 
1057 void Wcs::shiftReferencePixel(double dx, double dy) {
1058  assert(_wcsInfo);
1059  _wcsInfo->crpix[0] += dx;
1060  _wcsInfo->crpix[1] += dy;
1061 
1062  // tells libwcs to invalidate cached data, since transformation has been modified
1063  _wcsInfo->flag = 0;
1064 }
1065 
1066 /*
1067  * Now WCSA, pixel coordinates, but allowing for X0 and Y0
1068  */
1069 namespace detail {
1070 
1075  std::string const& wcsName,
1076  int const x0,
1077  int const y0
1078  ) {
1080 
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");
1089 
1090  return wcsMetaData;
1091 }
1102  int x0 = 0; // Our value of X0
1103  int y0 = 0; // Our value of Y0
1104 
1105  try {
1106  //
1107  // Only use WCS if CRPIX[12] == 1 and CRVAL[12] is present
1108  //
1109  if (metadata->getAsDouble("CRPIX1" + wcsName) == 1 &&
1110  metadata->getAsDouble("CRPIX2" + wcsName) == 1) {
1111  x0 = metadata->getAsInt("CRVAL1" + wcsName);
1112  y0 = metadata->getAsInt("CRVAL2" + wcsName);
1113  //
1114  // OK, we've got it. Remove it from the header
1115  //
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);
1124  }
1125  } catch (pex::exceptions::NotFoundError&) {
1126  ; // OK, not present
1127  }
1128 
1129  return geom::Point2I(x0, y0);
1130 }
1131 
1145  std::shared_ptr<daf::base::PropertySet> wcsMetadata = wcs->getFitsMetadata();
1146  std::vector<std::string> paramNames = wcsMetadata->paramNames();
1147  paramNames.push_back("CDELT1");
1148  paramNames.push_back("CDELT2");
1149  paramNames.push_back("LTV1");
1150  paramNames.push_back("LTV2");
1151  paramNames.push_back("PC001001");
1152  paramNames.push_back("PC001002");
1153  paramNames.push_back("PC002001");
1154  paramNames.push_back("PC002002");
1155  for (std::vector<std::string>::const_iterator ptr = paramNames.begin(); ptr != paramNames.end(); ++ptr) {
1156  metadata->remove(*ptr);
1157  }
1158 
1159  return 0; // would be ncard if remove returned a status
1160 }
1161 }
1162 
1163 // -------------------------------------------------------------------------------------------------
1164 //
1165 // XYTransformFromWcsPair
1166 
1168  : XYTransform(), _dst(dst), _src(src), _isSameSkySystem(dst->isSameSkySystem(*src)) {}
1169 
1171  return std::make_shared<XYTransformFromWcsPair>(_dst->clone(), _src->clone());
1172 }
1173 
1175  if (_isSameSkySystem) {
1176  // high performance branch; no coordinate conversion required
1177  afw::geom::Angle sky0, sky1;
1178  _src->pixelToSky(pixel[0], pixel[1], sky0, sky1);
1179  return _dst->skyToPixel(sky0, sky1);
1180  } else {
1181  std::shared_ptr<afw::coord::Coord const> coordPtr{_src->pixelToSky(pixel)};
1182  return _dst->skyToPixel(*coordPtr);
1183  }
1184 }
1185 
1187  if (_isSameSkySystem) {
1188  // high performance branch; no coordinate conversion required
1189  afw::geom::Angle sky0, sky1;
1190  _dst->pixelToSky(pixel[0], pixel[1], sky0, sky1);
1191  return _src->skyToPixel(sky0, sky1);
1192  } else {
1193  std::shared_ptr<afw::coord::Coord const> coordPtr{_dst->pixelToSky(pixel)};
1194  return _src->skyToPixel(*coordPtr);
1195  }
1196 }
1197 
1199  // just swap src, dst
1200  return std::make_shared<XYTransformFromWcsPair>(_src, _dst);
1201 }
1202 }
1203 }
1204 } // end image
double getAsDouble(std::string const &name) const
Get the last value for any arithmetic property name (possibly hierarchical).
Definition: PropertySet.cc:288
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.
Definition: Coord.cc:846
virtual ~Wcs()
Definition: Wcs.cc:525
void _setWcslibParams()
Definition: Wcs.cc:68
bool operator==(Wcs const &other) const
Definition: Wcs.cc:468
static std::shared_ptr< lsst::daf::base::PropertyList > generatePropertySet(lsst::afw::image::Wcs const &wcs)
virtual bool isPersistable() const
Whether the Wcs is persistable using afw::table::io archives.
Definition: Wcs.cc:1011
bool _isIcrs() const
Definition: Wcs.h:371
geom::Angle pixelScale() const
Returns the pixel scale [Angle/pixel].
Definition: Wcs.cc:697
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...
Definition: Wcs.cc:741
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
Definition: ApCorrMap.cc:71
afw::coord::CoordSystem getCoordSystem() const
Definition: Wcs.h:229
WcsFactory(std::string const &name)
Definition: Wcs.cc:922
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...
Definition: Wcs.cc:910
virtual std::shared_ptr< afw::geom::XYTransform > invert() const
returns a "deep inverse" in this sense that the forward+inverse transforms do not share state default...
Definition: Wcs.cc:1198
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).
Definition: PropertySet.cc:554
T front(T... args)
geom::Point2D skyToPixel(geom::Angle sky1, geom::Angle sky2) const
Convert from sky coordinates (e.g.
Definition: Wcs.cc:745
A custom container class for records, based on std::vector.
Definition: Catalog.h:97
Class for storing ordered metadata with comments.
Definition: PropertyList.h:73
AngleUnit constexpr radians
constant with units of radians
Definition: Angle.h:87
bool _mayBePersistable() const
Perform basic checks on whether *this might be persistable.
Definition: Wcs.cc:1003
A base class for factory classes used to reconstruct objects from records.
Definition: Persistable.h:217
std::string formatFitsProperties(lsst::daf::base::PropertySet const &prop, std::set< std::string > const &excludeNames={})
Format a PropertySet into a single FITS header string.
Definition: Utils.cc:267
int getAsInt(std::string const &name) const
Get the last value for a bool/char/short/int property name (possibly hierarchical).
Definition: PropertySet.cc:231
int stripWcsKeywords(std::shared_ptr< lsst::daf::base::PropertySet > const &metadata, std::shared_ptr< Wcs const > const &wcs)
Definition: Wcs.cc:1143
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:902
T calloc(T... args)
T endl(T... args)
AngleUnit constexpr degrees
constant with units of degrees
Definition: Angle.h:88
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:868
virtual Point2D forwardTransform(Point2D const &pixel) const
virtuals for forward and reverse transforms
Definition: Wcs.cc:1174
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...
Definition: Wcs.cc:644
lsst::afw::geom::Point2D GeomPoint
Definition: Wcs.cc:57
table::PointKey< double > crval
Definition: Wcs.cc:934
bool isSameSkySystem(Wcs const &wcs) const
Return true if a WCS has the same coordinate system and equinox as this one.
Definition: Wcs.cc:796
STL namespace.
tbl::Key< int > wcs
geom::Point2I getImageXY0FromMetadata(daf::base::PropertySet &metadata, std::string const &wcsName, bool strip=false)
table::Key< std::string > radesys
Definition: Wcs.cc:940
T end(T... args)
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:543
Implementation of the WCS standard for a any projection.
Definition: Wcs.h:104
friend class WcsFactory
Definition: Wcs.h:359
lsst::daf::base::PropertyList PropertyList
Definition: Wcs.cc:56
table::Key< std::string > ctype2
Definition: Wcs.cc:938
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition: BaseRecord.h:123
A class used to convert scalar POD types such as double to Angle.
Definition: Angle.h:52
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.
Definition: Wcs.cc:666
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:346
std::map< Citizen const *, CitizenInfo > table
Definition: Citizen.h:90
STL class.
double constexpr PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:21
A 2D linear coordinate transformation.
lsst::daf::base::PropertySet PropertySet
Definition: Wcs.cc:55
bool isFlipped() const
Does the Wcs follow the convention of North=Up, East=Left?
Definition: Wcs.cc:648
T push_back(T... args)
afw::table::PointKey< int > dimensions
Definition: GaussianPsf.cc:44
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.
Definition: Wcs.cc:160
A class representing an angle.
Definition: Angle.h:102
T strncpy(T... args)
Wcs & operator=(const Wcs &)
A base class for image defects.
Definition: cameraGeom.dox:3
std::shared_ptr< coord::Coord > pixelToSky(double pix1, double pix2) const
Convert from pixel position to sky coordinates (e.g.
Definition: Wcs.cc:827
virtual bool _isSubset(Wcs const &other) const
Definition: Wcs.cc:504
T strlen(T... args)
geom::LinearTransform getLinearTransform() const
Return the linear part of the Wcs, the CD matrix in FITS-speak, as an AffineTransform.
Definition: Wcs.cc:916
Point< double, 2 > Point2D
Definition: Point.h:300
std::shared_ptr< lsst::daf::base::PropertyList > createTrivialWcsAsPropertySet(std::string const &wcsName, int const x0=0, int const y0=0)
Definition: Wcs.cc:1074
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.
Definition: Point.h:122
An affine coordinate transformation consisting of a linear transformation and an offset.
table::Key< double > equinox
Definition: Wcs.cc:939
T strcmp(T... args)
const int fitsToLsstPixels
Definition: Wcs.cc:78
#define CHECK_NULLS(a, b)
Definition: Wcs.cc:493
virtual void shiftReferencePixel(double dx, double dy)
Move the pixel reference position by (dx, dy)
Definition: Wcs.cc:1057
const int STRLEN
Definition: Wcs.cc:61
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:134
Wcs()
Construct an invalid Wcs given no arguments.
Definition: Wcs.cc:84
table::Key< int > b
Include files required for standard LSST Exception handling.
int countFitsHeaderCards(lsst::daf::base::PropertySet const &prop)
Definition: Utils.cc:279
Schema getSchema() const
Return the Schema that holds this record&#39;s fields and keys.
Definition: BaseRecord.h:52
lsst::afw::geom::Point2D getPixelOrigin() const
Returns CRPIX (corrected to LSST convention).
Definition: Wcs.cc:542
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...
Definition: Wcs.cc:846
void _initWcs()
Set some internal variables that we need to refer to.
Definition: Wcs.cc:114
double x
geom::Point2D skyToIntermediateWorldCoord(coord::Coord const &coord) const
Convert from sky coordinates (e.g.
Definition: Wcs.cc:747
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
std::shared_ptr< lsst::afw::coord::Coord > getSkyOrigin() const
Returns CRVAL. This need not be the centre of the image.
Definition: Wcs.cc:537
virtual Point2D reverseTransform(Point2D const &pixel) const
Definition: Wcs.cc:1186
virtual void rotateImageBy90(int nQuarter, lsst::afw::geom::Extent2I dimensions) const
Rotate image by nQuarter times 90 degrees.
Definition: Wcs.cc:591
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
table::Schema schema
Definition: Wcs.cc:933
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:982
std::shared_ptr< Coord > convert(CoordSystem system, double epoch=2000) const
Convert to a specified Coord type at a specified epoch.
Definition: Coord.cc:471
T size(T... args)
XYTransformFromWcsPair(std::shared_ptr< Wcs const > dst, std::shared_ptr< Wcs const > src)
Definition: Wcs.cc:1167
bool _skyAxesSwapped
if true then the sky axes are swapped
Definition: Wcs.h:453
#define LSST_EXCEPT(type,...)
Create an exception with a given type and message and optionally other arguments (dependent on the ty...
Definition: Exception.h:46
A vector of catalogs used by Persistable.
Definition: CatalogVector.h:29
T malloc(T... args)
Base class for all records.
Definition: BaseRecord.h:29
int _wcshdrCtrl
Controls messages to stderr from wcshdr (0 for none); see wcshdr.h for details.
Definition: Wcs.h:450
int _relax
Degree of permissiveness for wcspih (0 for strict); see wcshdr.h for details.
Definition: Wcs.h:448
Citizen(const std::type_info &)
Definition: Citizen.cc:174
T begin(T... args)
virtual void pixelToSkyImpl(double pixel1, double pixel2, geom::Angle skyTmp[2]) const
Worker routine for pixelToSky.
Definition: Wcs.cc:803
virtual std::shared_ptr< afw::geom::XYTransform > clone() const
The following methods are needed to devirtualize the XYTransform parent class.
Definition: Wcs.cc:1170
std::shared_ptr< lsst::afw::coord::Coord > CoordPtr
Definition: Wcs.cc:58
Class for storing generic metadata.
Definition: PropertySet.h:73
T c_str(T... args)
std::shared_ptr< Wcs const > _dst
Definition: Wcs.h:521
virtual geom::Point2D skyToPixelImpl(geom::Angle sky1, geom::Angle sky2) const
Worker routine for skyToPixel.
Definition: Wcs.cc:699
int _wcsfixCtrl
Do potentially unsafe translations of non-standard unit strings? 0/1 = no/yes.
Definition: Wcs.h:449
T substr(T... args)
virtual void write(OutputArchiveHandle &handle) const
Write the object to one or more catalogs.
Definition: Wcs.cc:984
table::Key< std::string > cunit1
Definition: Wcs.cc:941
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:112
table::Key< int > a
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...
Definition: Wcs.cc:875
table::Key< table::Array< double > > cd
Definition: Wcs.cc:936
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.
Definition: Wcs.cc:1041
double getEquinox() const
Definition: Wcs.cc:789
const int lsstToFitsPixels
Definition: Wcs.cc:77
T quiet_NaN(T... args)
T sqrt(T... args)
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
Point< int, 2 > Point2I
Definition: Point.h:297
Eigen::Matrix2d getCDMatrix() const
Returns the CD matrix.
Definition: Wcs.cc:550
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:341
A multi-catalog archive object used to load table::io::Persistable objects.
Definition: InputArchive.h:31
constexpr double asDegrees() const noexcept
Return an Angle&#39;s value in degrees.
Definition: Angle.h:144
This is the base class for spherical coordinates.
Definition: Coord.h:63
const double PixelZeroPos
position of center of pixel 0
Definition: ImageUtils.h:44
ImageT val
Definition: CR.cc:158
virtual std::shared_ptr< Wcs > clone(void) const
Definition: Wcs.cc:531
T isspace(T... args)
AffineTransform const invert() const
Return the inverse transform.
table::Key< std::string > ctype1
Definition: Wcs.cc:937
std::shared_ptr< Wcs const > _src
Definition: Wcs.h:522
struct wcsprm * _wcsInfo
Definition: Wcs.h:446
coord::CoordSystem _coordSystem
Definition: Wcs.h:452
virtual void flipImage(int flipLR, int flipTB, lsst::afw::geom::Extent2I dimensions) const
Flip CD matrix around the y-axis.
Definition: Wcs.cc:568
An exception thrown when an InputArchive&#39;s contents do not make sense.
Definition: Persistable.h:39
table::Key< std::string > cunit2
Definition: Wcs.cc:942
std::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:472
virtual std::string getPersistenceName() const
Return the unique name used to persist this object and look up its factory.
Definition: Wcs.cc:980
table::PointKey< double > crpix
Definition: Wcs.cc:935