LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
SkyWcs.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2017 AURA/LSST.
4  *
5  * This product includes software developed by the
6  * LSST Project (http://www.lsst.org/).
7  *
8  * This program is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the LSST License Statement and
19  * the GNU General Public License along with this program. If not,
20  * see <http://www.lsstcorp.org/LegalNotices/>.
21  */
22 
23 #include <cmath>
24 #include <cstdint>
25 #include <exception>
26 #include <memory>
27 #include <vector>
28 
29 #include "astshim.h"
30 
31 #include "lsst/geom/Angle.h"
32 #include "lsst/geom/Point.h"
33 #include "lsst/geom/SpherePoint.h"
38 #include "lsst/afw/geom/wcsUtils.h"
39 #include "lsst/afw/geom/SkyWcs.h"
41 #include "lsst/pex/exceptions.h"
43 
44 namespace lsst {
45 namespace afw {
46 
49 
50 namespace geom {
51 namespace {
52 
53 int const SERIALIZATION_VERSION = 1;
54 
55 // TIGHT_FITS_TOL is used by getFitsMetadata to determine if a WCS can accurately be represented as a FITS
56 // WCS. It specifies the maximum departure from linearity (in pixels) allowed on either axis of the mapping
57 // from pixel coordinates to Intermediate World Coordinates over a range of 100 x 100 pixels
58 // (or the region specified by NAXIS[12], if provided, but we do not pass this to AST as of 2018-01-17).
59 // For more information,
60 // see FitsTol in the AST manual http://starlink.eao.hawaii.edu/devdocs/sun211.htx/sun211.html
61 double const TIGHT_FITS_TOL = 0.0001;
62 
63 class SkyWcsPersistenceHelper {
64 public:
65  table::Schema schema;
66  table::Key<table::Array<std::uint8_t>> wcs;
67 
68  static SkyWcsPersistenceHelper const& get() {
69  static SkyWcsPersistenceHelper instance;
70  return instance;
71  }
72 
73  // No copying
74  SkyWcsPersistenceHelper(const SkyWcsPersistenceHelper&) = delete;
75  SkyWcsPersistenceHelper& operator=(const SkyWcsPersistenceHelper&) = delete;
76 
77  // No moving
78  SkyWcsPersistenceHelper(SkyWcsPersistenceHelper&&) = delete;
79  SkyWcsPersistenceHelper& operator=(SkyWcsPersistenceHelper&&) = delete;
80 
81 private:
82  SkyWcsPersistenceHelper()
83  : schema(),
84  wcs(schema.addField<table::Array<std::uint8_t>>("wcs", "wcs string representation", "")) {}
85 };
86 
87 class SkyWcsFactory : public table::io::PersistableFactory {
88 public:
89  explicit SkyWcsFactory(std::string const& name) : table::io::PersistableFactory(name) {}
90 
91  std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
92  CatalogVector const& catalogs) const override {
93  SkyWcsPersistenceHelper const& keys = SkyWcsPersistenceHelper::get();
94  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
95  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
96  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
97  table::BaseRecord const& record = catalogs.front().front();
98  std::string stringRep = formatters::bytesToString(record.get(keys.wcs));
99  return SkyWcs::readString(stringRep);
100  }
101 };
102 
103 std::string getSkyWcsPersistenceName() { return "SkyWcs"; }
104 
105 SkyWcsFactory registration(getSkyWcsPersistenceName());
106 
107 ast::FrameDict makeSkyWcsFrameDict(TransformPoint2ToPoint2 const& pixelsToFieldAngle,
108  lsst::geom::Angle const& orientation, bool flipX,
110  std::string const& projection = "TAN") {
111  auto const orientationAndFlipXMatrix = makeCdMatrix(1 * lsst::geom::degrees, orientation, flipX);
112  auto const initialWcs =
113  makeSkyWcs(lsst::geom::Point2D(0, 0), crval, orientationAndFlipXMatrix, projection);
114  auto const initialFrameDict = initialWcs->getFrameDict();
115  auto const iwcToSkyMap = initialFrameDict->getMapping("IWC", "SKY");
116  auto const pixelFrame = initialFrameDict->getFrame("PIXELS");
117  auto const iwcFrame = initialFrameDict->getFrame("IWC");
118  auto const skyFrame = initialFrameDict->getFrame("SKY");
119  // Field angle is in radians and is aligned to focal plane x and y;
120  // IWC is basically the same thing, but in degrees and with rotation and flipX applied
121  ndarray::Array<double, 2, 2> fieldAngleToIwcNdArray = ndarray::allocate(2, 2);
122  asEigenMatrix(fieldAngleToIwcNdArray) = orientationAndFlipXMatrix * 180.0 / lsst::geom::PI;
123  auto const pixelsToFieldAngleMap = pixelsToFieldAngle.getMapping();
124  auto const fieldAngleToIwcMap = ast::MatrixMap(fieldAngleToIwcNdArray);
125  auto const pixelsToIwcMap = pixelsToFieldAngleMap->then(fieldAngleToIwcMap);
126  auto finalFrameDict = ast::FrameDict(*pixelFrame, pixelsToIwcMap, *iwcFrame);
127  finalFrameDict.addFrame("IWC", *iwcToSkyMap, *skyFrame);
128  return finalFrameDict;
129 }
130 
131 } // namespace
132 
134  bool flipX) {
135  Eigen::Matrix2d cdMatrix;
136  double orientRad = orientation.asRadians();
137  double scaleDeg = scale.asDegrees();
138  double xmult = flipX ? 1.0 : -1.0;
139  cdMatrix(0, 0) = std::cos(orientRad) * scaleDeg * xmult;
140  cdMatrix(0, 1) = std::sin(orientRad) * scaleDeg;
141  cdMatrix(1, 0) = -std::sin(orientRad) * scaleDeg * xmult;
142  cdMatrix(1, 1) = std::cos(orientRad) * scaleDeg;
143  return cdMatrix;
144 }
145 
147  auto const dstInverse = dst.getTransform()->inverted();
148  return src.getTransform()->then(*dstInverse);
149 }
150 
151 SkyWcs::SkyWcs(daf::base::PropertySet& metadata, bool strip)
152  : SkyWcs(detail::readLsstSkyWcs(metadata, strip)) {}
153 
154 SkyWcs::SkyWcs(ast::FrameDict const& frameDict) : SkyWcs(_checkFrameDict(frameDict)) {}
155 
156 bool SkyWcs::operator==(SkyWcs const& other) const { return writeString() == other.writeString(); }
157 
159  // Compute pixVec containing the pixel position and two nearby points
160  // (use a vector so all three points can be converted to sky in a single call)
161  double const side = 1.0;
163  pixel,
164  pixel + lsst::geom::Extent2D(side, 0),
165  pixel + lsst::geom::Extent2D(0, side),
166  };
167 
168  auto skyVec = pixelToSky(pixVec);
169 
170  // Work in 3-space to avoid RA wrapping and pole issues
171  auto skyLL = skyVec[0].getVector();
172  auto skyDx = skyVec[1].getVector() - skyLL;
173  auto skyDy = skyVec[2].getVector() - skyLL;
174 
175  // Compute pixel scale in radians = sqrt(pixel area in radians^2)
176  // pixel area in radians^2 = area of parallelogram with sides skyDx, skyDy = |skyDx cross skyDy|
177  // Use squared norm to avoid two square roots
178  double skyAreaSq = skyDx.cross(skyDy).getSquaredNorm();
179  return (std::pow(skyAreaSq, 0.25) / side) * lsst::geom::radians;
180 }
181 
183  // CRVAL is stored as the SkyRef property of the sky frame (the current frame of the SkyWcs)
184  auto skyFrame = std::dynamic_pointer_cast<ast::SkyFrame>(
185  getFrameDict()->getFrame(ast::FrameDict::CURRENT, false)); // false: do not copy
186  if (!skyFrame) {
187  throw LSST_EXCEPT(pex::exceptions::LogicError, "Current frame is not a SkyFrame");
188  }
189  auto const crvalRad = skyFrame->getSkyRef();
190  return lsst::geom::SpherePoint(crvalRad[0] * lsst::geom::radians, crvalRad[1] * lsst::geom::radians);
191 }
192 
193 Eigen::Matrix2d SkyWcs::getCdMatrix(lsst::geom::Point2D const& pixel) const {
194  auto const pixelToIwc = getFrameDict()->getMapping(ast::FrameSet::BASE, "IWC");
195  auto const pixelToIwcTransform = TransformPoint2ToPoint2(*pixelToIwc);
196  return pixelToIwcTransform.getJacobian(pixel);
197 }
198 
199 Eigen::Matrix2d SkyWcs::getCdMatrix() const { return getCdMatrix(getPixelOrigin()); }
200 
202  auto const crval = pixelToSky(pixel);
203  auto const cdMatrix = getCdMatrix(pixel);
204  auto metadata = makeSimpleWcsMetadata(pixel, crval, cdMatrix);
205  return std::make_shared<SkyWcs>(*metadata);
206 }
207 
209  auto newToOldPixel = TransformPoint2ToPoint2(ast::ShiftMap({-shift[0], -shift[1]}));
210  return makeModifiedWcs(newToOldPixel, *this, true);
211 }
212 
214  // Make a FrameSet that maps from GRID to SKY; GRID = the base frame (PIXELS or ACTUAL_PIXELS) + 1
215  auto const gridToPixel = ast::ShiftMap({-1.0, -1.0});
216  auto thisDict = getFrameDict();
217  auto const pixelToIwc = thisDict->getMapping(ast::FrameSet::BASE, "IWC");
218  auto const iwcToSky = thisDict->getMapping("IWC", "SKY");
219  auto const gridToSky = gridToPixel.then(*pixelToIwc).then(*iwcToSky);
220  ast::FrameSet frameSet(ast::Frame(2, "Domain=GRID"), gridToSky, *thisDict->getFrame("SKY", false));
221 
222  // Write frameSet to a FitsChan and extract the metadata
224  os << "Encoding=FITS-WCS, CDMatrix=1, FitsAxisOrder=<copy>, FitsTol=" << TIGHT_FITS_TOL;
225  ast::StringStream strStream;
226  ast::FitsChan fitsChan(strStream, os.str());
227  int const nObjectsWritten = fitsChan.write(frameSet);
229  if (nObjectsWritten == 0) {
230  if (precise) {
232  "Could not represent this SkyWcs using FITS-WCS metadata");
233  } else {
234  // An exact representation could not be written, so try to write a local TAN approximation;
235  // set precise true to avoid an infinite loop, should something go wrong
236  auto tanWcs = getTanWcs(getPixelOrigin());
237  header = tanWcs->getFitsMetadata(true);
238  }
239  } else {
240  header = detail::getPropertyListFromFitsChan(fitsChan);
241  }
242 
243  // Remove DATE-OBS, MJD-OBS: AST writes these if the EQUINOX is set, but we set them via other mechanisms.
244  header->remove("DATE-OBS");
245  header->remove("MJD-OBS");
246 
247  // If CD matrix is present, explicitly set any missing entries to zero, as a convenience to the user
248  bool const hasCd11 = header->exists("CD1_1");
249  bool const hasCd12 = header->exists("CD1_2");
250  bool const hasCd21 = header->exists("CD2_1");
251  bool const hasCd22 = header->exists("CD2_2");
252  if (hasCd11 || hasCd12 || hasCd21 || hasCd22) {
253  if (!hasCd11) header->set("CD1_1", 0.0, "Transformation matrix element");
254  if (!hasCd12) header->set("CD1_2", 0.0, "Transformation matrix element");
255  if (!hasCd21) header->set("CD2_1", 0.0, "Transformation matrix element");
256  if (!hasCd22) header->set("CD2_2", 0.0, "Transformation matrix element");
257  }
258 
259  return header;
260 }
261 
263 
264 bool SkyWcs::isFits() const {
265  try {
266  getFitsMetadata(true);
267  } catch (const lsst::pex::exceptions::RuntimeError&) {
268  return false;
269  } catch (const std::runtime_error&) {
270  return false;
271  }
272  return true;
273 }
274 
276  lsst::geom::AngleUnit const& skyUnit) const {
277  return _linearizePixelToSky(skyToPixel(coord), coord, skyUnit);
278 }
280  lsst::geom::AngleUnit const& skyUnit) const {
281  return _linearizePixelToSky(pix, pixelToSky(pix), skyUnit);
282 }
283 
285  lsst::geom::AngleUnit const& skyUnit) const {
286  return _linearizeSkyToPixel(skyToPixel(coord), coord, skyUnit);
287 }
288 
290  lsst::geom::AngleUnit const& skyUnit) const {
291  return _linearizeSkyToPixel(pix, pixelToSky(pix), skyUnit);
292 }
293 
294 std::string SkyWcs::getShortClassName() { return "SkyWcs"; };
295 
296 bool SkyWcs::isFlipped() const {
297  double det = getCdMatrix().determinant();
298  if (det == 0) {
299  throw(LSST_EXCEPT(pex::exceptions::RuntimeError, "CD matrix is singular"));
300  }
301  return (det > 0);
302 }
303 
305  int version;
306  is >> version;
307  if (version != 1) {
308  throw LSST_EXCEPT(pex::exceptions::TypeError, "Unsupported version " + std::to_string(version));
309  }
310  std::string shortClassName;
311  is >> shortClassName;
312  if (shortClassName != SkyWcs::getShortClassName()) {
314  os << "Class name in stream " << shortClassName << " != " << SkyWcs::getShortClassName();
316  }
318  is >> hasFitsApproximation;
319  auto astStream = ast::Stream(&is, nullptr);
320  auto astObjectPtr = ast::Channel(astStream).read();
321  auto frameSet = std::dynamic_pointer_cast<ast::FrameSet>(astObjectPtr);
322  if (!frameSet) {
324  os << "The AST serialization was read as a " << astObjectPtr->getClassName()
325  << " instead of a FrameSet";
327  }
328  ast::FrameDict frameDict(*frameSet);
329  return std::make_shared<SkyWcs>(frameDict);
330 }
331 
333  std::istringstream is(str);
334  return SkyWcs::readStream(is);
335 }
336 
338  os << SERIALIZATION_VERSION << " " << SkyWcs::getShortClassName() << " " << hasFitsApproximation() << " ";
339  getFrameDict()->show(os, false); // false = do not write comments
340 }
341 
344  writeStream(os);
345  return os.str();
346 }
347 
349  return std::make_unique<SkyWcs>(*this);
350 }
351 
354  if (isFits()) {
355  os << "FITS standard SkyWcs:";
356  } else {
357  os << "Non-standard SkyWcs (Frames: ";
358  // Print the frames in index order (frames are numbered from 1).
359  std::string delimiter = "";
360  for (size_t i = 1; i <= getFrameDict()->getAllDomains().size(); ++i) {
361  os << delimiter << getFrameDict()->getFrame(i)->getDomain();
362  delimiter = ", ";
363  }
364  os << "): ";
365  }
366  std::string delimiter = "\n";
367  os << delimiter << "Sky Origin: " << getSkyOrigin();
368  os << delimiter << "Pixel Origin: " << getPixelOrigin();
369  os << delimiter << "Pixel Scale: " << getPixelScale().asArcseconds() << " arcsec/pixel";
370  return os.str();
371 }
372 
373 bool SkyWcs::equals(typehandling::Storable const& other) const noexcept {
374  return singleClassEquals(*this, other);
375 }
376 
377 std::string SkyWcs::getPersistenceName() const { return getSkyWcsPersistenceName(); }
378 
379 std::string SkyWcs::getPythonModule() const { return "lsst.afw.geom"; }
380 
381 void SkyWcs::write(OutputArchiveHandle& handle) const {
382  SkyWcsPersistenceHelper const& keys = SkyWcsPersistenceHelper::get();
383  table::BaseCatalog cat = handle.makeCatalog(keys.schema);
385  record->set(keys.wcs, formatters::stringToBytes(writeString()));
386  handle.saveCatalog(cat);
387 }
388 
390  : _frameDict(frameDict), _transform(), _pixelOrigin(), _pixelScaleAtOrigin(0 * lsst::geom::radians) {
391  _computeCache();
392 };
393 
394 std::shared_ptr<ast::FrameDict> SkyWcs::_checkFrameDict(ast::FrameDict const& frameDict) const {
395  // Check that each frame is present and has the right type and number of axes
396  std::vector<std::string> const domainNames = {"ACTUAL_PIXELS", "PIXELS", "IWC", "SKY"};
397  for (auto const& domainName : domainNames) {
398  if (frameDict.hasDomain(domainName)) {
399  auto const frame = frameDict.getFrame(domainName, false);
400  if (frame->getNAxes() != 2) {
402  "Frame " + domainName + " has " + std::to_string(frame->getNAxes()) +
403  " axes instead of 2");
404  }
405  auto desiredClassName = domainName == "SKY" ? "SkyFrame" : "Frame";
406  if (frame->getClassName() != desiredClassName) {
408  "Frame " + domainName + " is of type " + frame->getClassName() +
409  " instead of " + desiredClassName);
410  }
411  } else if (domainName != "ACTUAL_PIXELS") {
412  // This is a required frame
414  "No frame with domain " + domainName + " found");
415  }
416  }
417 
418  // The base frame must have domain "PIXELS" or "ACTUAL_PIXELS"
419  auto baseDomain = frameDict.getFrame(ast::FrameSet::BASE, false)->getDomain();
420  if (baseDomain != "ACTUAL_PIXELS" && baseDomain != "PIXELS") {
422  "Base frame has domain " + baseDomain + " instead of PIXELS or ACTUAL_PIXELS");
423  }
424 
425  // The current frame must have domain "SKY"
426  auto currentDomain = frameDict.getFrame(ast::FrameSet::CURRENT, false)->getDomain();
427  if (currentDomain != "SKY") {
429  "Current frame has domain " + currentDomain + " instead of SKY");
430  }
431 
432  return frameDict.copy();
433 }
434 
435 lsst::geom::AffineTransform SkyWcs::_linearizePixelToSky(lsst::geom::Point2D const& pix00,
436  lsst::geom::SpherePoint const& coord,
437  lsst::geom::AngleUnit const& skyUnit) const {
438  // Figure out the (0, 0), (0, 1), and (1, 0) ra/dec coordinates of the corners
439  // of a square drawn in pixel. It'd be better to center the square at sky00,
440  // but that would involve another conversion between sky and pixel coordinates
441  const double side = 1.0; // length of the square's sides in pixels
442  auto const sky00 = coord.getPosition(skyUnit);
443  auto const dsky10 = coord.getTangentPlaneOffset(pixelToSky(pix00 + lsst::geom::Extent2D(side, 0)));
444  auto const dsky01 = coord.getTangentPlaneOffset(pixelToSky(pix00 + lsst::geom::Extent2D(0, side)));
445 
446  Eigen::Matrix2d m;
447  m(0, 0) = dsky10.first.asAngularUnits(skyUnit) / side;
448  m(0, 1) = dsky01.first.asAngularUnits(skyUnit) / side;
449  m(1, 0) = dsky10.second.asAngularUnits(skyUnit) / side;
450  m(1, 1) = dsky01.second.asAngularUnits(skyUnit) / side;
451 
452  Eigen::Vector2d sky00v;
453  sky00v << sky00.getX(), sky00.getY();
454  Eigen::Vector2d pix00v;
455  pix00v << pix00.getX(), pix00.getY();
456  // return lsst::geom::AffineTransform(m, lsst::geom::Extent2D(sky00v - m * pix00v));
457  return lsst::geom::AffineTransform(m, (sky00v - m * pix00v));
458 }
459 
460 lsst::geom::AffineTransform SkyWcs::_linearizeSkyToPixel(lsst::geom::Point2D const& pix00,
461  lsst::geom::SpherePoint const& coord,
462  lsst::geom::AngleUnit const& skyUnit) const {
463  lsst::geom::AffineTransform inverse = _linearizePixelToSky(pix00, coord, skyUnit);
464  return inverse.inverted();
465 }
466 
467 std::shared_ptr<SkyWcs> makeFlippedWcs(SkyWcs const& wcs, bool flipLR, bool flipTB,
468  lsst::geom::Point2D const& center) {
469  double const dx = 1000; // any "reasonable" number of pixels will do
470  std::vector<double> inLL = {center[0] - dx, center[1] - dx};
471  std::vector<double> inUR = {center[0] + dx, center[1] + dx};
472  std::vector<double> outLL(inLL);
473  std::vector<double> outUR(inUR);
474  if (flipLR) {
475  outLL[0] = inUR[0];
476  outUR[0] = inLL[0];
477  }
478  if (flipTB) {
479  outLL[1] = inUR[1];
480  outUR[1] = inLL[1];
481  }
482  auto const flipPix = TransformPoint2ToPoint2(ast::WinMap(inLL, inUR, outLL, outUR));
483  return makeModifiedWcs(flipPix, wcs, true);
484 }
485 
487  bool modifyActualPixels) {
488  auto const pixelMapping = pixelTransform.getMapping();
489  auto oldFrameDict = wcs.getFrameDict();
490  bool const hasActualPixels = oldFrameDict->hasDomain("ACTUAL_PIXELS");
491  auto const pixelFrame = oldFrameDict->getFrame("PIXELS", false);
492  auto const iwcFrame = oldFrameDict->getFrame("IWC", false);
493  auto const skyFrame = oldFrameDict->getFrame("SKY", false);
494  auto const oldPixelToIwc = oldFrameDict->getMapping("PIXELS", "IWC");
495  auto const iwcToSky = oldFrameDict->getMapping("IWC", "SKY");
496 
497  std::shared_ptr<ast::FrameDict> newFrameDict;
498  std::shared_ptr<ast::Mapping> newPixelToIwc;
499  if (hasActualPixels) {
500  auto const actualPixelFrame = oldFrameDict->getFrame("ACTUAL_PIXELS", false);
501  auto const oldActualPixelToPixels = oldFrameDict->getMapping("ACTUAL_PIXELS", "PIXELS");
502  std::shared_ptr<ast::Mapping> newActualPixelsToPixels;
503  if (modifyActualPixels) {
504  newActualPixelsToPixels = pixelMapping->then(*oldActualPixelToPixels).simplified();
505  newPixelToIwc = oldPixelToIwc;
506  } else {
507  newActualPixelsToPixels = oldActualPixelToPixels;
508  newPixelToIwc = pixelMapping->then(*oldPixelToIwc).simplified();
509  }
510  newFrameDict =
511  std::make_shared<ast::FrameDict>(*actualPixelFrame, *newActualPixelsToPixels, *pixelFrame);
512  newFrameDict->addFrame("PIXELS", *newPixelToIwc, *iwcFrame);
513  } else {
514  newPixelToIwc = pixelMapping->then(*oldPixelToIwc).simplified();
515  newFrameDict = std::make_shared<ast::FrameDict>(*pixelFrame, *newPixelToIwc, *iwcFrame);
516  }
517  newFrameDict->addFrame("IWC", *iwcToSky, *skyFrame);
518  return std::make_shared<SkyWcs>(*newFrameDict);
519 }
520 
522  return std::make_shared<SkyWcs>(metadata, strip);
523 }
524 
526  Eigen::Matrix2d const& cdMatrix, std::string const& projection) {
527  auto metadata = makeSimpleWcsMetadata(crpix, crval, cdMatrix, projection);
528  return std::make_shared<SkyWcs>(*metadata);
529 }
530 
532  lsst::geom::Angle const& orientation, bool flipX,
533  lsst::geom::SpherePoint const& boresight, std::string const& projection) {
534  auto frameDict = makeSkyWcsFrameDict(pixelsToFieldAngle, orientation, flipX, boresight, projection);
535  return std::make_shared<SkyWcs>(frameDict);
536 }
537 
539  Eigen::Matrix2d const& cdMatrix, Eigen::MatrixXd const& sipA,
540  Eigen::MatrixXd const& sipB) {
541  auto metadata = makeTanSipMetadata(crpix, crval, cdMatrix, sipA, sipB);
542  return std::make_shared<SkyWcs>(*metadata);
543 }
544 
546  Eigen::Matrix2d const& cdMatrix, Eigen::MatrixXd const& sipA,
547  Eigen::MatrixXd const& sipB, Eigen::MatrixXd const& sipAp,
548  Eigen::MatrixXd const& sipBp) {
549  auto metadata = makeTanSipMetadata(crpix, crval, cdMatrix, sipA, sipB, sipAp, sipBp);
550  return std::make_shared<SkyWcs>(*metadata);
551 }
552 
554  bool simplify) {
555  auto iwcToSky = wcs.getFrameDict()->getMapping("IWC", "SKY");
556  return std::make_shared<TransformPoint2ToSpherePoint>(*iwcToSky, simplify);
557 }
558 
560  auto pixelToIwc = wcs.getFrameDict()->getMapping(ast::FrameSet::BASE, "IWC");
561  return std::make_shared<TransformPoint2ToPoint2>(*pixelToIwc, simplify);
562 }
563 
565  os << wcs.toString();
566  return os;
567 };
568 
569 } // namespace geom
570 } // namespace afw
571 } // namespace lsst
table::Key< std::string > name
Definition: Amplifier.cc:116
table::PointKey< int > pixel
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
std::shared_ptr< RecordT > src
Definition: Match.cc:48
table::PointKey< double > crpix
Definition: OldWcs.cc:129
table::PointKey< double > crval
Definition: OldWcs.cc:128
std::ostream * os
Definition: Schema.cc:557
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:66
table::Schema schema
Definition: SkyWcs.cc:65
int m
Definition: SpanSet.cc:48
table::Key< table::Array< std::uint8_t > > frameSet
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
Channel provides input/output of AST objects.
Definition: Channel.h:60
std::shared_ptr< Object > read()
Read an object from a channel.
Definition: Channel.cc:52
int write(Object const &object)
Write an object to a channel.
Definition: Channel.cc:61
A specialized form of Channel which reads and writes FITS header cards.
Definition: FitsChan.h:202
A FrameSet whose frames can be referenced by domain name.
Definition: FrameDict.h:67
bool hasDomain(std::string const &domain) const
Return True if a frame in this FrameDict has the specified domain.
Definition: FrameDict.h:201
std::shared_ptr< FrameDict > copy() const
Return a deep copy of this object.
Definition: FrameDict.h:123
std::set< std::string > getAllDomains() const
Get the domain names for all contained Frames (excluding frames with empty or defaulted domain names)...
Definition: FrameDict.cc:45
std::shared_ptr< Mapping > getMapping(int from, std::string const &to) const
Variant of FrameSet::getMapping with the second frame specified by domain.
Definition: FrameDict.h:162
std::shared_ptr< Frame > getFrame(std::string const &domain, bool copy=true) const
Variant of getFrame(int, bool) where the frame is specified by domain name.
Definition: FrameDict.h:151
Frame is used to represent a coordinate system.
Definition: Frame.h:157
A FrameSet consists of a set of one or more Frames (which describe coordinate systems),...
Definition: FrameSet.h:99
static constexpr int CURRENT
index of current frame
Definition: FrameSet.h:105
static constexpr int BASE
index of base frame
Definition: FrameSet.h:104
SeriesMap then(Mapping const &next) const
Return a series compound mapping this(first(input)) containing shallow copies of the original.
Definition: Mapping.cc:37
std::shared_ptr< Mapping > simplified() const
Return a simplied version of the mapping (which may be a compound Mapping such as a CmpMap).
Definition: Mapping.h:248
MatrixMap is a form of Mapping which performs a general linear transformation.
Definition: MatrixMap.h:42
void show(std::ostream &os, bool showComments=true) const
Print a textual description the object to an ostream.
Definition: Object.cc:158
ShiftMap is a linear Mapping which shifts each axis by a specified constant value.
Definition: ShiftMap.h:40
A stream for ast::Channel.
Definition: Stream.h:41
String-based source and sink for channels.
Definition: Stream.h:180
A WinMap is a linear Mapping which transforms a rectangular window in one coordinate system into a si...
Definition: WinMap.h:45
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
Eigen::Matrix2d getCdMatrix() const
Get the 2x2 CD matrix at the pixel origin.
Definition: SkyWcs.cc:199
std::string writeString() const
Serialize this SkyWcs to a string, using the same format as writeStream.
Definition: SkyWcs.cc:342
lsst::geom::SpherePoint pixelToSky(lsst::geom::Point2D const &pixel) const
Compute sky position(s) from pixel position(s)
Definition: SkyWcs.h:334
std::shared_ptr< SkyWcs > copyAtShiftedPixelOrigin(lsst::geom::Extent2D const &shift) const
Return a copy of this SkyWcs with the pixel origin shifted by the specified amount.
Definition: SkyWcs.cc:208
static std::shared_ptr< SkyWcs > readStream(std::istream &is)
Deserialize a SkyWcs from an input stream.
Definition: SkyWcs.cc:304
std::shared_ptr< SkyWcs > getTanWcs(lsst::geom::Point2D const &pixel) const
Get a local TAN WCS approximation to this WCS at the specified pixel position.
Definition: SkyWcs.cc:201
lsst::geom::Angle getPixelScale() const
Get the pixel scale at the pixel origin.
Definition: SkyWcs.h:208
static std::shared_ptr< SkyWcs > readString(std::string &str)
Deserialize a SkyWcs from a string, using the same format as readStream.
Definition: SkyWcs.cc:332
lsst::geom::SpherePoint getSkyOrigin() const
Get the sky origin, the celestial fiducial point.
Definition: SkyWcs.cc:182
std::string toString() const override
Create a string representation of this object.
Definition: SkyWcs.cc:352
lsst::geom::Point2D getPixelOrigin() const
Get the pixel origin, in pixels, using the LSST convention.
Definition: SkyWcs.h:215
lsst::geom::AffineTransform linearizePixelToSky(lsst::geom::SpherePoint const &coord, lsst::geom::AngleUnit const &skyUnit) const
Return the local linear approximation to pixelToSky at a point given in sky coordinates.
Definition: SkyWcs.cc:275
bool operator==(SkyWcs const &other) const
Equality is based on the string representations being equal.
Definition: SkyWcs.cc:156
lsst::geom::AffineTransform linearizeSkyToPixel(lsst::geom::SpherePoint const &coord, lsst::geom::AngleUnit const &skyUnit) const
Return the local linear approximation to skyToPixel at a point given in sky coordinates.
Definition: SkyWcs.cc:284
SkyWcs(SkyWcs const &)=default
lsst::geom::Point2D skyToPixel(lsst::geom::SpherePoint const &sky) const
Compute pixel position(s) from sky position(s)
Definition: SkyWcs.h:349
std::string getPythonModule() const override
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
Definition: SkyWcs.cc:379
bool isFlipped() const
Does the WCS follow the convention of North=Up, East=Left?
Definition: SkyWcs.cc:296
std::shared_ptr< daf::base::PropertyList > getFitsMetadata(bool precise=false) const
Return the WCS as FITS WCS metadata.
Definition: SkyWcs.cc:213
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
Definition: SkyWcs.cc:377
bool hasFitsApproximation() const
Does this SkyWcs have an approximate SkyWcs that can be represented as standard FITS WCS?
Definition: SkyWcs.h:364
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
Definition: SkyWcs.cc:381
bool equals(typehandling::Storable const &other) const noexcept override
Compare this object to another Storable.
Definition: SkyWcs.cc:373
std::shared_ptr< const TransformPoint2ToSpherePoint > getTransform() const
Get a TransformPoint2ToSpherePoint that transforms pixels to sky in the forward direction and sky to ...
Definition: SkyWcs.h:257
void writeStream(std::ostream &os) const
Serialize this SkyWcs to an output stream.
Definition: SkyWcs.cc:337
std::shared_ptr< typehandling::Storable > cloneStorable() const override
Create a new SkyWcs that is a copy of this one.
Definition: SkyWcs.cc:348
static std::string getShortClassName()
Definition: SkyWcs.cc:294
std::shared_ptr< const ast::FrameDict > getFrameDict() const
Get the contained FrameDict.
Definition: SkyWcs.cc:262
bool isFits() const
Return true getFitsMetadata(true) will succeed, false if not.
Definition: SkyWcs.cc:264
Transform LSST spatial data, such as lsst::geom::Point2D and lsst::geom::SpherePoint,...
Definition: Transform.h:68
std::shared_ptr< const ast::Mapping > getMapping() const
Get the contained mapping.
Definition: Transform.h:135
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:490
An object passed to Persistable::write to allow it to persist itself.
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
static std::shared_ptr< T > dynamicCast(std::shared_ptr< Persistable > const &ptr)
Dynamically cast a shared_ptr.
Definition: Persistable.cc:18
Interface supporting iteration over heterogenous containers.
Definition: Storable.h:58
Class for storing generic metadata.
Definition: PropertySet.h:66
An affine coordinate transformation consisting of a linear transformation and an offset.
AffineTransform const inverted() const
Return the inverse transform.
A class representing an angle.
Definition: Angle.h:127
constexpr double asArcseconds() const noexcept
Return an Angle's value in arcseconds.
Definition: Angle.h:178
A class used to convert scalar POD types such as double to Angle.
Definition: Angle.h:70
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
std::pair< Angle, Angle > getTangentPlaneOffset(SpherePoint const &other) const
Get the offset from a tangent plane centered at this point to another point.
Definition: SpherePoint.cc:214
Point2D getPosition(AngleUnit unit) const noexcept
Return longitude, latitude as a Point2D object.
Definition: SpherePoint.cc:129
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
Reports errors from accepting an object of an unexpected or inappropriate type.
Definition: Runtime.h:167
T cos(T... args)
bool strip
Definition: fits.cc:911
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:108
ndarray::Array< std::uint8_t, 1, 1 > stringToBytes(std::string const &str)
Encode a std::string as a vector of uint8.
Definition: Utils.cc:161
std::string bytesToString(ndarray::Array< std::uint8_t const, 1, 1 > const &bytes)
Decode a std::string from a vector of uint8 returned by stringToBytes.
Definition: Utils.cc:172
std::shared_ptr< ast::FrameDict > readLsstSkyWcs(daf::base::PropertySet &metadata, bool strip=true)
Read an LSST celestial WCS FrameDict from a FITS header.
std::shared_ptr< daf::base::PropertyList > getPropertyListFromFitsChan(ast::FitsChan &fitsChan)
Copy values from an AST FitsChan into a PropertyList.
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition: SkyWcs.cc:521
std::shared_ptr< daf::base::PropertyList > makeSimpleWcsMetadata(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval, Eigen::Matrix2d const &cdMatrix, std::string const &projection="TAN")
Make FITS metadata for a simple FITS WCS (one with no distortion).
Definition: wcsUtils.cc:169
std::shared_ptr< TransformPoint2ToPoint2 > getPixelToIntermediateWorldCoords(SkyWcs const &wcs, bool simplify=true)
Return a transform from pixel coordinates to intermediate world coordinates.
Definition: SkyWcs.cc:559
std::shared_ptr< SkyWcs > makeTanSipWcs(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval, Eigen::Matrix2d const &cdMatrix, Eigen::MatrixXd const &sipA, Eigen::MatrixXd const &sipB)
Construct a TAN-SIP SkyWcs with forward SIP distortion terms and an iterative inverse.
Definition: SkyWcs.cc:538
std::ostream & operator<<(std::ostream &os, GenericEndpoint const &endpoint)
Print "GenericEndpoint(_n_)" to the ostream where _n_ is the number of axes, e.g. "GenericAxes(4)".
Definition: Endpoint.cc:239
Eigen::Matrix2d makeCdMatrix(lsst::geom::Angle const &scale, lsst::geom::Angle const &orientation=0 *lsst::geom::degrees, bool flipX=false)
Make a WCS CD matrix.
Definition: SkyWcs.cc:133
std::shared_ptr< TransformPoint2ToSpherePoint > getIntermediateWorldCoordsToSky(SkyWcs const &wcs, bool simplify=true)
Return a transform from intermediate world coordinates to sky.
Definition: SkyWcs.cc:553
std::shared_ptr< daf::base::PropertyList > makeTanSipMetadata(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval, Eigen::Matrix2d const &cdMatrix, Eigen::MatrixXd const &sipA, Eigen::MatrixXd const &sipB)
Make metadata for a TAN-SIP WCS without inverse matrices.
Definition: wcsUtils.cc:194
std::shared_ptr< SkyWcs > makeModifiedWcs(TransformPoint2ToPoint2 const &pixelTransform, SkyWcs const &wcs, bool modifyActualPixels)
Create a new SkyWcs whose pixels are transformed by pixelTransform, as described below.
Definition: SkyWcs.cc:486
std::shared_ptr< TransformPoint2ToPoint2 > makeWcsPairTransform(SkyWcs const &src, SkyWcs const &dst)
A Transform obtained by putting two SkyWcs objects "back to back".
Definition: SkyWcs.cc:146
Transform< Point2Endpoint, Point2Endpoint > TransformPoint2ToPoint2
Definition: Transform.h:300
std::shared_ptr< SkyWcs > makeFlippedWcs(SkyWcs const &wcs, bool flipLR, bool flipTB, lsst::geom::Point2D const &center)
Return a copy of a FITS-WCS with pixel positions flipped around a specified center.
Definition: SkyWcs.cc:467
FilterProperty & operator=(FilterProperty const &)=default
lsst::geom::SpherePoint SpherePoint
Definition: misc.h:34
constexpr double PI
The ratio of a circle's circumference to diameter.
Definition: Angle.h:39
Extent< double, 2 > Extent2D
Definition: Extent.h:400
constexpr AngleUnit degrees
constant with units of degrees
Definition: Angle.h:109
constexpr AngleUnit radians
constant with units of radians
Definition: Angle.h:108
int orientation(UnitVector3d const &a, UnitVector3d const &b, UnitVector3d const &c)
orientation computes and returns the orientations of 3 unit vectors a, b and c.
Definition: orientation.cc:135
A base class for image defects.
STL namespace.
T pow(T... args)
T sin(T... args)
T size(T... args)
T to_string(T... args)
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override
Definition: warpExposure.cc:0