LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
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 <algorithm>
24 #include <cmath>
25 #include <cstdint>
26 #include <exception>
27 #include <memory>
28 #include <ostream>
29 #include <sstream>
30 #include <vector>
31 
32 #include "astshim.h"
33 
34 #include "lsst/geom/Angle.h"
35 #include "lsst/geom/Point.h"
36 #include "lsst/geom/SpherePoint.h"
38 #include "lsst/afw/table.h"
43 #include "lsst/afw/geom/wcsUtils.h"
44 #include "lsst/afw/geom/SkyWcs.h"
46 #include "lsst/pex/exceptions.h"
48 
49 namespace lsst {
50 namespace afw {
51 
54 
55 namespace geom {
56 namespace {
57 
58 int const SERIALIZATION_VERSION = 1;
59 
60 // TIGHT_FITS_TOL is used by getFitsMetadata to determine if a WCS can accurately be represented as a FITS
61 // WCS. It specifies the maximum departure from linearity (in pixels) allowed on either axis of the mapping
62 // from pixel coordinates to Intermediate World Coordinates over a range of 100 x 100 pixels
63 // (or the region specified by NAXIS[12], if provided, but we do not pass this to AST as of 2018-01-17).
64 // For more information,
65 // see FitsTol in the AST manual http://starlink.eao.hawaii.edu/devdocs/sun211.htx/sun211.html
66 double const TIGHT_FITS_TOL = 0.0001;
67 
68 class SkyWcsPersistenceHelper {
69 public:
70  table::Schema schema;
71  table::Key<table::Array<std::uint8_t>> wcs;
72 
73  static SkyWcsPersistenceHelper const& get() {
74  static SkyWcsPersistenceHelper instance;
75  return instance;
76  }
77 
78  // No copying
79  SkyWcsPersistenceHelper(const SkyWcsPersistenceHelper&) = delete;
80  SkyWcsPersistenceHelper& operator=(const SkyWcsPersistenceHelper&) = delete;
81 
82  // No moving
83  SkyWcsPersistenceHelper(SkyWcsPersistenceHelper&&) = delete;
84  SkyWcsPersistenceHelper& operator=(SkyWcsPersistenceHelper&&) = delete;
85 
86 private:
87  SkyWcsPersistenceHelper()
88  : schema(),
89  wcs(schema.addField<table::Array<std::uint8_t>>("wcs", "wcs string representation", "")) {}
90 };
91 
92 class SkyWcsFactory : public table::io::PersistableFactory {
93 public:
94  explicit SkyWcsFactory(std::string const& name) : table::io::PersistableFactory(name) {}
95 
96  std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
97  CatalogVector const& catalogs) const override {
98  SkyWcsPersistenceHelper const& keys = SkyWcsPersistenceHelper::get();
99  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
100  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
101  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
102  table::BaseRecord const& record = catalogs.front().front();
103  std::string stringRep = formatters::bytesToString(record.get(keys.wcs));
104  return SkyWcs::readString(stringRep);
105  }
106 };
107 
108 std::string getSkyWcsPersistenceName() { return "SkyWcs"; }
109 
110 SkyWcsFactory registration(getSkyWcsPersistenceName());
111 
112 ast::FrameDict makeSkyWcsFrameDict(TransformPoint2ToPoint2 const& pixelsToFieldAngle,
113  lsst::geom::Angle const& orientation, bool flipX,
115  std::string const& projection = "TAN") {
116  auto const orientationAndFlipXMatrix = makeCdMatrix(1 * lsst::geom::degrees, orientation, flipX);
117  auto const initialWcs =
118  makeSkyWcs(lsst::geom::Point2D(0, 0), crval, orientationAndFlipXMatrix, projection);
119  auto const initialFrameDict = initialWcs->getFrameDict();
120  auto const iwcToSkyMap = initialFrameDict->getMapping("IWC", "SKY");
121  auto const pixelFrame = initialFrameDict->getFrame("PIXELS");
122  auto const iwcFrame = initialFrameDict->getFrame("IWC");
123  auto const skyFrame = initialFrameDict->getFrame("SKY");
124  // Field angle is in radians and is aligned to focal plane x and y;
125  // IWC is basically the same thing, but in degrees and with rotation and flipX applied
126  ndarray::Array<double, 2, 2> fieldAngleToIwcNdArray = ndarray::allocate(2, 2);
127  asEigenMatrix(fieldAngleToIwcNdArray) = orientationAndFlipXMatrix * 180.0 / lsst::geom::PI;
128  auto const pixelsToFieldAngleMap = pixelsToFieldAngle.getMapping();
129  auto const fieldAngleToIwcMap = ast::MatrixMap(fieldAngleToIwcNdArray);
130  auto const pixelsToIwcMap = pixelsToFieldAngleMap->then(fieldAngleToIwcMap);
131  auto finalFrameDict = ast::FrameDict(*pixelFrame, pixelsToIwcMap, *iwcFrame);
132  finalFrameDict.addFrame("IWC", *iwcToSkyMap, *skyFrame);
133  return finalFrameDict;
134 }
135 
136 } // namespace
137 
139  bool flipX) {
140  Eigen::Matrix2d cdMatrix;
141  double orientRad = orientation.asRadians();
142  double scaleDeg = scale.asDegrees();
143  double xmult = flipX ? 1.0 : -1.0;
144  cdMatrix(0, 0) = std::cos(orientRad) * scaleDeg * xmult;
145  cdMatrix(0, 1) = std::sin(orientRad) * scaleDeg;
146  cdMatrix(1, 0) = -std::sin(orientRad) * scaleDeg * xmult;
147  cdMatrix(1, 1) = std::cos(orientRad) * scaleDeg;
148  return cdMatrix;
149 }
150 
152  auto const dstInverse = dst.getTransform()->inverted();
153  return src.getTransform()->then(*dstInverse);
154 }
155 
156 SkyWcs::SkyWcs(daf::base::PropertySet& metadata, bool strip)
157  : SkyWcs(detail::readLsstSkyWcs(metadata, strip)) {}
158 
159 SkyWcs::SkyWcs(ast::FrameDict const& frameDict) : SkyWcs(_checkFrameDict(frameDict)) {}
160 
161 bool SkyWcs::operator==(SkyWcs const& other) const { return writeString() == other.writeString(); }
162 
164  // Compute pixVec containing the pixel position and two nearby points
165  // (use a vector so all three points can be converted to sky in a single call)
166  double const side = 1.0;
168  pixel,
169  pixel + lsst::geom::Extent2D(side, 0),
170  pixel + lsst::geom::Extent2D(0, side),
171  };
172 
173  auto skyVec = pixelToSky(pixVec);
174 
175  // Work in 3-space to avoid RA wrapping and pole issues
176  auto skyLL = skyVec[0].getVector();
177  auto skyDx = skyVec[1].getVector() - skyLL;
178  auto skyDy = skyVec[2].getVector() - skyLL;
179 
180  // Compute pixel scale in radians = sqrt(pixel area in radians^2)
181  // pixel area in radians^2 = area of parallelogram with sides skyDx, skyDy = |skyDx cross skyDy|
182  // Use squared norm to avoid two square roots
183  double skyAreaSq = skyDx.cross(skyDy).getSquaredNorm();
184  return (std::pow(skyAreaSq, 0.25) / side) * lsst::geom::radians;
185 }
186 
188  // CRVAL is stored as the SkyRef property of the sky frame (the current frame of the SkyWcs)
189  auto skyFrame = std::dynamic_pointer_cast<ast::SkyFrame>(
190  getFrameDict()->getFrame(ast::FrameDict::CURRENT, false)); // false: do not copy
191  if (!skyFrame) {
192  throw LSST_EXCEPT(pex::exceptions::LogicError, "Current frame is not a SkyFrame");
193  }
194  auto const crvalRad = skyFrame->getSkyRef();
195  return lsst::geom::SpherePoint(crvalRad[0] * lsst::geom::radians, crvalRad[1] * lsst::geom::radians);
196 }
197 
198 Eigen::Matrix2d SkyWcs::getCdMatrix(lsst::geom::Point2D const& pixel) const {
199  auto const pixelToIwc = getFrameDict()->getMapping(ast::FrameSet::BASE, "IWC");
200  auto const pixelToIwcTransform = TransformPoint2ToPoint2(*pixelToIwc);
201  return pixelToIwcTransform.getJacobian(pixel);
202 }
203 
204 Eigen::Matrix2d SkyWcs::getCdMatrix() const { return getCdMatrix(getPixelOrigin()); }
205 
207  auto const crval = pixelToSky(pixel);
208  auto const cdMatrix = getCdMatrix(pixel);
209  auto metadata = makeSimpleWcsMetadata(pixel, crval, cdMatrix);
210  return std::make_shared<SkyWcs>(*metadata);
211 }
212 
214  auto newToOldPixel = TransformPoint2ToPoint2(ast::ShiftMap({-shift[0], -shift[1]}));
215  return makeModifiedWcs(newToOldPixel, *this, true);
216 }
217 
219  // Make a FrameSet that maps from GRID to SKY; GRID = the base frame (PIXELS or ACTUAL_PIXELS) + 1
220  auto const gridToPixel = ast::ShiftMap({-1.0, -1.0});
221  auto thisDict = getFrameDict();
222  auto const pixelToIwc = thisDict->getMapping(ast::FrameSet::BASE, "IWC");
223  auto const iwcToSky = thisDict->getMapping("IWC", "SKY");
224  auto const gridToSky = gridToPixel.then(*pixelToIwc).then(*iwcToSky);
225  ast::FrameSet frameSet(ast::Frame(2, "Domain=GRID"), gridToSky, *thisDict->getFrame("SKY", false));
226 
227  // Write frameSet to a FitsChan and extract the metadata
229  os << "Encoding=FITS-WCS, CDMatrix=1, FitsAxisOrder=<copy>, FitsTol=" << TIGHT_FITS_TOL;
230  ast::StringStream strStream;
231  ast::FitsChan fitsChan(strStream, os.str());
232  int const nObjectsWritten = fitsChan.write(frameSet);
234  if (nObjectsWritten == 0) {
235  if (precise) {
237  "Could not represent this SkyWcs using FITS-WCS metadata");
238  } else {
239  // An exact representation could not be written, so try to write a local TAN approximation;
240  // set precise true to avoid an infinite loop, should something go wrong
241  auto tanWcs = getTanWcs(getPixelOrigin());
242  header = tanWcs->getFitsMetadata(true);
243  }
244  } else {
245  header = detail::getPropertyListFromFitsChan(fitsChan);
246  }
247 
248  // Remove DATE-OBS, MJD-OBS: AST writes these if the EQUINOX is set, but we set them via other mechanisms.
249  header->remove("DATE-OBS");
250  header->remove("MJD-OBS");
251 
252  // If CD matrix is present, explicitly set any missing entries to zero, as a convenience to the user
253  bool const hasCd11 = header->exists("CD1_1");
254  bool const hasCd12 = header->exists("CD1_2");
255  bool const hasCd21 = header->exists("CD2_1");
256  bool const hasCd22 = header->exists("CD2_2");
257  if (hasCd11 || hasCd12 || hasCd21 || hasCd22) {
258  if (!hasCd11) header->set("CD1_1", 0.0, "Transformation matrix element");
259  if (!hasCd12) header->set("CD1_2", 0.0, "Transformation matrix element");
260  if (!hasCd21) header->set("CD2_1", 0.0, "Transformation matrix element");
261  if (!hasCd22) header->set("CD2_2", 0.0, "Transformation matrix element");
262  }
263 
264  return header;
265 }
266 
268 
269 bool SkyWcs::isFits() const {
270  try {
271  getFitsMetadata(true);
272  } catch (const lsst::pex::exceptions::RuntimeError&) {
273  return false;
274  } catch (const std::runtime_error&) {
275  return false;
276  }
277  return true;
278 }
279 
281  lsst::geom::AngleUnit const& skyUnit) const {
282  return _linearizePixelToSky(skyToPixel(coord), coord, skyUnit);
283 }
285  lsst::geom::AngleUnit const& skyUnit) const {
286  return _linearizePixelToSky(pix, pixelToSky(pix), skyUnit);
287 }
288 
290  lsst::geom::AngleUnit const& skyUnit) const {
291  return _linearizeSkyToPixel(skyToPixel(coord), coord, skyUnit);
292 }
293 
295  lsst::geom::AngleUnit const& skyUnit) const {
296  return _linearizeSkyToPixel(pix, pixelToSky(pix), skyUnit);
297 }
298 
299 std::string SkyWcs::getShortClassName() { return "SkyWcs"; };
300 
301 bool SkyWcs::isFlipped() const {
302  double det = getCdMatrix().determinant();
303  if (det == 0) {
304  throw(LSST_EXCEPT(pex::exceptions::RuntimeError, "CD matrix is singular"));
305  }
306  return (det > 0);
307 }
308 
310  int version;
311  is >> version;
312  if (version != 1) {
313  throw LSST_EXCEPT(pex::exceptions::TypeError, "Unsupported version " + std::to_string(version));
314  }
315  std::string shortClassName;
316  is >> shortClassName;
317  if (shortClassName != SkyWcs::getShortClassName()) {
319  os << "Class name in stream " << shortClassName << " != " << SkyWcs::getShortClassName();
321  }
323  is >> hasFitsApproximation;
324  auto astStream = ast::Stream(&is, nullptr);
325  auto astObjectPtr = ast::Channel(astStream).read();
326  auto frameSet = std::dynamic_pointer_cast<ast::FrameSet>(astObjectPtr);
327  if (!frameSet) {
329  os << "The AST serialization was read as a " << astObjectPtr->getClassName()
330  << " instead of a FrameSet";
332  }
333  ast::FrameDict frameDict(*frameSet);
334  return std::make_shared<SkyWcs>(frameDict);
335 }
336 
338  std::istringstream is(str);
339  return SkyWcs::readStream(is);
340 }
341 
343  os << SERIALIZATION_VERSION << " " << SkyWcs::getShortClassName() << " " << hasFitsApproximation() << " ";
344  getFrameDict()->show(os, false); // false = do not write comments
345 }
346 
349  writeStream(os);
350  return os.str();
351 }
352 
354  return std::make_unique<SkyWcs>(*this);
355 }
356 
359  if (isFits()) {
360  os << "FITS standard SkyWcs:";
361  } else {
362  os << "Non-standard SkyWcs (Frames: ";
363  // Print the frames in index order (frames are numbered from 1).
364  std::string delimiter = "";
365  for (size_t i = 1; i <= getFrameDict()->getAllDomains().size(); ++i) {
366  os << delimiter << getFrameDict()->getFrame(i)->getDomain();
367  delimiter = ", ";
368  }
369  os << "): ";
370  }
371  std::string delimiter = "\n";
372  os << delimiter << "Sky Origin: " << getSkyOrigin();
373  os << delimiter << "Pixel Origin: " << getPixelOrigin();
374  os << delimiter << "Pixel Scale: " << getPixelScale().asArcseconds() << " arcsec/pixel";
375  return os.str();
376 }
377 
378 bool SkyWcs::equals(typehandling::Storable const& other) const noexcept {
379  return singleClassEquals(*this, other);
380 }
381 
382 std::string SkyWcs::getPersistenceName() const { return getSkyWcsPersistenceName(); }
383 
384 std::string SkyWcs::getPythonModule() const { return "lsst.afw.geom"; }
385 
386 void SkyWcs::write(OutputArchiveHandle& handle) const {
387  SkyWcsPersistenceHelper const& keys = SkyWcsPersistenceHelper::get();
388  table::BaseCatalog cat = handle.makeCatalog(keys.schema);
390  record->set(keys.wcs, formatters::stringToBytes(writeString()));
391  handle.saveCatalog(cat);
392 }
393 
395  : _frameDict(frameDict), _transform(), _pixelOrigin(), _pixelScaleAtOrigin(0 * lsst::geom::radians) {
396  _computeCache();
397 };
398 
399 std::shared_ptr<ast::FrameDict> SkyWcs::_checkFrameDict(ast::FrameDict const& frameDict) const {
400  // Check that each frame is present and has the right type and number of axes
401  std::vector<std::string> const domainNames = {"ACTUAL_PIXELS", "PIXELS", "IWC", "SKY"};
402  for (auto const& domainName : domainNames) {
403  if (frameDict.hasDomain(domainName)) {
404  auto const frame = frameDict.getFrame(domainName, false);
405  if (frame->getNAxes() != 2) {
407  "Frame " + domainName + " has " + std::to_string(frame->getNAxes()) +
408  " axes instead of 2");
409  }
410  auto desiredClassName = domainName == "SKY" ? "SkyFrame" : "Frame";
411  if (frame->getClassName() != desiredClassName) {
413  "Frame " + domainName + " is of type " + frame->getClassName() +
414  " instead of " + desiredClassName);
415  }
416  } else if (domainName != "ACTUAL_PIXELS") {
417  // This is a required frame
419  "No frame with domain " + domainName + " found");
420  }
421  }
422 
423  // The base frame must have domain "PIXELS" or "ACTUAL_PIXELS"
424  auto baseDomain = frameDict.getFrame(ast::FrameSet::BASE, false)->getDomain();
425  if (baseDomain != "ACTUAL_PIXELS" && baseDomain != "PIXELS") {
427  "Base frame has domain " + baseDomain + " instead of PIXELS or ACTUAL_PIXELS");
428  }
429 
430  // The current frame must have domain "SKY"
431  auto currentDomain = frameDict.getFrame(ast::FrameSet::CURRENT, false)->getDomain();
432  if (currentDomain != "SKY") {
434  "Current frame has domain " + currentDomain + " instead of SKY");
435  }
436 
437  return frameDict.copy();
438 }
439 
440 lsst::geom::AffineTransform SkyWcs::_linearizePixelToSky(lsst::geom::Point2D const& pix00,
441  lsst::geom::SpherePoint const& coord,
442  lsst::geom::AngleUnit const& skyUnit) const {
443  // Figure out the (0, 0), (0, 1), and (1, 0) ra/dec coordinates of the corners
444  // of a square drawn in pixel. It'd be better to center the square at sky00,
445  // but that would involve another conversion between sky and pixel coordinates
446  const double side = 1.0; // length of the square's sides in pixels
447  auto const sky00 = coord.getPosition(skyUnit);
448  auto const dsky10 = coord.getTangentPlaneOffset(pixelToSky(pix00 + lsst::geom::Extent2D(side, 0)));
449  auto const dsky01 = coord.getTangentPlaneOffset(pixelToSky(pix00 + lsst::geom::Extent2D(0, side)));
450 
451  Eigen::Matrix2d m;
452  m(0, 0) = dsky10.first.asAngularUnits(skyUnit) / side;
453  m(0, 1) = dsky01.first.asAngularUnits(skyUnit) / side;
454  m(1, 0) = dsky10.second.asAngularUnits(skyUnit) / side;
455  m(1, 1) = dsky01.second.asAngularUnits(skyUnit) / side;
456 
457  Eigen::Vector2d sky00v;
458  sky00v << sky00.getX(), sky00.getY();
459  Eigen::Vector2d pix00v;
460  pix00v << pix00.getX(), pix00.getY();
461  // return lsst::geom::AffineTransform(m, lsst::geom::Extent2D(sky00v - m * pix00v));
462  return lsst::geom::AffineTransform(m, (sky00v - m * pix00v));
463 }
464 
465 lsst::geom::AffineTransform SkyWcs::_linearizeSkyToPixel(lsst::geom::Point2D const& pix00,
466  lsst::geom::SpherePoint const& coord,
467  lsst::geom::AngleUnit const& skyUnit) const {
468  lsst::geom::AffineTransform inverse = _linearizePixelToSky(pix00, coord, skyUnit);
469  return inverse.inverted();
470 }
471 
472 std::shared_ptr<SkyWcs> makeFlippedWcs(SkyWcs const& wcs, bool flipLR, bool flipTB,
473  lsst::geom::Point2D const& center) {
474  double const dx = 1000; // any "reasonable" number of pixels will do
475  std::vector<double> inLL = {center[0] - dx, center[1] - dx};
476  std::vector<double> inUR = {center[0] + dx, center[1] + dx};
477  std::vector<double> outLL(inLL);
478  std::vector<double> outUR(inUR);
479  if (flipLR) {
480  outLL[0] = inUR[0];
481  outUR[0] = inLL[0];
482  }
483  if (flipTB) {
484  outLL[1] = inUR[1];
485  outUR[1] = inLL[1];
486  }
487  auto const flipPix = TransformPoint2ToPoint2(ast::WinMap(inLL, inUR, outLL, outUR));
488  return makeModifiedWcs(flipPix, wcs, true);
489 }
490 
492  bool modifyActualPixels) {
493  auto const pixelMapping = pixelTransform.getMapping();
494  auto oldFrameDict = wcs.getFrameDict();
495  bool const hasActualPixels = oldFrameDict->hasDomain("ACTUAL_PIXELS");
496  auto const pixelFrame = oldFrameDict->getFrame("PIXELS", false);
497  auto const iwcFrame = oldFrameDict->getFrame("IWC", false);
498  auto const skyFrame = oldFrameDict->getFrame("SKY", false);
499  auto const oldPixelToIwc = oldFrameDict->getMapping("PIXELS", "IWC");
500  auto const iwcToSky = oldFrameDict->getMapping("IWC", "SKY");
501 
502  std::shared_ptr<ast::FrameDict> newFrameDict;
503  std::shared_ptr<ast::Mapping> newPixelToIwc;
504  if (hasActualPixels) {
505  auto const actualPixelFrame = oldFrameDict->getFrame("ACTUAL_PIXELS", false);
506  auto const oldActualPixelToPixels = oldFrameDict->getMapping("ACTUAL_PIXELS", "PIXELS");
507  std::shared_ptr<ast::Mapping> newActualPixelsToPixels;
508  if (modifyActualPixels) {
509  newActualPixelsToPixels = pixelMapping->then(*oldActualPixelToPixels).simplified();
510  newPixelToIwc = oldPixelToIwc;
511  } else {
512  newActualPixelsToPixels = oldActualPixelToPixels;
513  newPixelToIwc = pixelMapping->then(*oldPixelToIwc).simplified();
514  }
515  newFrameDict =
516  std::make_shared<ast::FrameDict>(*actualPixelFrame, *newActualPixelsToPixels, *pixelFrame);
517  newFrameDict->addFrame("PIXELS", *newPixelToIwc, *iwcFrame);
518  } else {
519  newPixelToIwc = pixelMapping->then(*oldPixelToIwc).simplified();
520  newFrameDict = std::make_shared<ast::FrameDict>(*pixelFrame, *newPixelToIwc, *iwcFrame);
521  }
522  newFrameDict->addFrame("IWC", *iwcToSky, *skyFrame);
523  return std::make_shared<SkyWcs>(*newFrameDict);
524 }
525 
527  return std::make_shared<SkyWcs>(metadata, strip);
528 }
529 
531  Eigen::Matrix2d const& cdMatrix, std::string const& projection) {
532  auto metadata = makeSimpleWcsMetadata(crpix, crval, cdMatrix, projection);
533  return std::make_shared<SkyWcs>(*metadata);
534 }
535 
537  lsst::geom::Angle const& orientation, bool flipX,
538  lsst::geom::SpherePoint const& boresight, std::string const& projection) {
539  auto frameDict = makeSkyWcsFrameDict(pixelsToFieldAngle, orientation, flipX, boresight, projection);
540  return std::make_shared<SkyWcs>(frameDict);
541 }
542 
544  Eigen::Matrix2d const& cdMatrix, Eigen::MatrixXd const& sipA,
545  Eigen::MatrixXd const& sipB) {
546  auto metadata = makeTanSipMetadata(crpix, crval, cdMatrix, sipA, sipB);
547  return std::make_shared<SkyWcs>(*metadata);
548 }
549 
551  Eigen::Matrix2d const& cdMatrix, Eigen::MatrixXd const& sipA,
552  Eigen::MatrixXd const& sipB, Eigen::MatrixXd const& sipAp,
553  Eigen::MatrixXd const& sipBp) {
554  auto metadata = makeTanSipMetadata(crpix, crval, cdMatrix, sipA, sipB, sipAp, sipBp);
555  return std::make_shared<SkyWcs>(*metadata);
556 }
557 
559  bool simplify) {
560  auto iwcToSky = wcs.getFrameDict()->getMapping("IWC", "SKY");
561  return std::make_shared<TransformPoint2ToSpherePoint>(*iwcToSky, simplify);
562 }
563 
565  auto pixelToIwc = wcs.getFrameDict()->getMapping(ast::FrameSet::BASE, "IWC");
566  return std::make_shared<TransformPoint2ToPoint2>(*pixelToIwc, simplify);
567 }
568 
570  os << wcs.toString();
571  return os;
572 };
573 
574 } // namespace geom
575 } // namespace afw
576 } // 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:131
table::PointKey< double > crval
Definition: OldWcs.cc:130
char delimiter
Definition: Schema.cc:193
ItemVariant const * other
Definition: Schema.cc:56
std::ostream * os
Definition: Schema.cc:746
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:71
table::Schema schema
Definition: SkyWcs.cc:70
int m
Definition: SpanSet.cc:49
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:204
std::string writeString() const
Serialize this SkyWcs to a string, using the same format as writeStream.
Definition: SkyWcs.cc:347
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:213
static std::shared_ptr< SkyWcs > readStream(std::istream &is)
Deserialize a SkyWcs from an input stream.
Definition: SkyWcs.cc:309
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:206
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:337
lsst::geom::SpherePoint getSkyOrigin() const
Get the sky origin, the celestial fiducial point.
Definition: SkyWcs.cc:187
std::string toString() const override
Create a string representation of this object.
Definition: SkyWcs.cc:357
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:280
bool operator==(SkyWcs const &other) const
Equality is based on the string representations being equal.
Definition: SkyWcs.cc:161
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:289
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:384
bool isFlipped() const
Does the WCS follow the convention of North=Up, East=Left?
Definition: SkyWcs.cc:301
std::shared_ptr< daf::base::PropertyList > getFitsMetadata(bool precise=false) const
Return the WCS as FITS WCS metadata.
Definition: SkyWcs.cc:218
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
Definition: SkyWcs.cc:382
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:386
bool equals(typehandling::Storable const &other) const noexcept override
Compare this object to another Storable.
Definition: SkyWcs.cc:378
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:342
std::shared_ptr< typehandling::Storable > cloneStorable() const override
Create a new SkyWcs that is a copy of this one.
Definition: SkyWcs.cc:353
static std::string getShortClassName()
Definition: SkyWcs.cc:299
std::shared_ptr< const ast::FrameDict > getFrameDict() const
Get the contained FrameDict.
Definition: SkyWcs.cc:267
bool isFits() const
Return true getFitsMetadata(true) will succeed, false if not.
Definition: SkyWcs.cc:269
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:485
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:67
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:162
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:173
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:526
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:172
std::shared_ptr< TransformPoint2ToPoint2 > getPixelToIntermediateWorldCoords(SkyWcs const &wcs, bool simplify=true)
Return a transform from pixel coordinates to intermediate world coordinates.
Definition: SkyWcs.cc:564
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:543
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:240
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:138
std::shared_ptr< TransformPoint2ToSpherePoint > getIntermediateWorldCoordsToSky(SkyWcs const &wcs, bool simplify=true)
Return a transform from intermediate world coordinates to sky.
Definition: SkyWcs.cc:558
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:197
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:491
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:151
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:472
FilterProperty & operator=(FilterProperty const &)=default
lsst::geom::SpherePoint SpherePoint
Definition: misc.h:35
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)