LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
LSSTDataManagementBasePackage
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)
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);
233  if (nObjectsWritten == 0) {
234  if (precise) {
236  "Could not represent this SkyWcs using FITS-WCS metadata");
237  } else {
238  // An exact representation could not be written, so try to write a local TAN approximation;
239  // set precise true to avoid an infinite loop, should something go wrong
240  auto tanWcs = getTanWcs(getPixelOrigin());
241  return tanWcs->getFitsMetadata(true);
242  }
243  }
245 }
246 
248 
249 bool SkyWcs::isFits() const {
250  try {
251  getFitsMetadata(true);
252  } catch (const lsst::pex::exceptions::RuntimeError&) {
253  return false;
254  } catch (const std::runtime_error&) {
255  return false;
256  }
257  return true;
258 }
259 
261  lsst::geom::AngleUnit const& skyUnit) const {
262  return _linearizePixelToSky(skyToPixel(coord), coord, skyUnit);
263 }
265  lsst::geom::AngleUnit const& skyUnit) const {
266  return _linearizePixelToSky(pix, pixelToSky(pix), skyUnit);
267 }
268 
270  lsst::geom::AngleUnit const& skyUnit) const {
271  return _linearizeSkyToPixel(skyToPixel(coord), coord, skyUnit);
272 }
273 
275  lsst::geom::AngleUnit const& skyUnit) const {
276  return _linearizeSkyToPixel(pix, pixelToSky(pix), skyUnit);
277 }
278 
279 std::string SkyWcs::getShortClassName() { return "SkyWcs"; };
280 
281 bool SkyWcs::isFlipped() const {
282  double det = getCdMatrix().determinant();
283  if (det == 0) {
284  throw(LSST_EXCEPT(pex::exceptions::RuntimeError, "CD matrix is singular"));
285  }
286  return (det > 0);
287 }
288 
290  int version;
291  is >> version;
292  if (version != 1) {
293  throw LSST_EXCEPT(pex::exceptions::TypeError, "Unsupported version " + std::to_string(version));
294  }
295  std::string shortClassName;
296  is >> shortClassName;
297  if (shortClassName != SkyWcs::getShortClassName()) {
299  os << "Class name in stream " << shortClassName << " != " << SkyWcs::getShortClassName();
301  }
303  is >> hasFitsApproximation;
304  auto astStream = ast::Stream(&is, nullptr);
305  auto astObjectPtr = ast::Channel(astStream).read();
306  auto frameSet = std::dynamic_pointer_cast<ast::FrameSet>(astObjectPtr);
307  if (!frameSet) {
309  os << "The AST serialization was read as a " << astObjectPtr->getClassName()
310  << " instead of a FrameSet";
312  }
314  return std::make_shared<SkyWcs>(frameDict);
315 }
316 
318  std::istringstream is(str);
319  return SkyWcs::readStream(is);
320 }
321 
323  os << SERIALIZATION_VERSION << " " << SkyWcs::getShortClassName() << " " << hasFitsApproximation() << " ";
324  getFrameDict()->show(os, false); // false = do not write comments
325 }
326 
329  writeStream(os);
330  return os.str();
331 }
332 
334  return std::make_unique<SkyWcs>(*this);
335 }
336 
339  if (isFits()) {
340  os << "FITS standard SkyWcs:";
341  } else {
342  os << "Non-standard SkyWcs (Frames: ";
343  // Print the frames in index order (frames are numbered from 1).
344  std::string delimiter = "";
345  for (size_t i = 1; i <= getFrameDict()->getAllDomains().size(); ++i) {
346  os << delimiter << getFrameDict()->getFrame(i)->getDomain();
347  delimiter = ", ";
348  }
349  os << "): ";
350  }
351  std::string delimiter = "\n";
352  os << delimiter << "Sky Origin: " << getSkyOrigin();
353  os << delimiter << "Pixel Origin: " << getPixelOrigin();
354  os << delimiter << "Pixel Scale: " << getPixelScale().asArcseconds() << " arcsec/pixel";
355  return os.str();
356 }
357 
358 bool SkyWcs::equals(typehandling::Storable const& other) const noexcept {
359  return singleClassEquals(*this, other);
360 }
361 
362 std::string SkyWcs::getPersistenceName() const { return getSkyWcsPersistenceName(); }
363 
364 std::string SkyWcs::getPythonModule() const { return "lsst.afw.geom"; }
365 
366 void SkyWcs::write(OutputArchiveHandle& handle) const {
367  SkyWcsPersistenceHelper const& keys = SkyWcsPersistenceHelper::get();
368  table::BaseCatalog cat = handle.makeCatalog(keys.schema);
370  record->set(keys.wcs, formatters::stringToBytes(writeString()));
371  handle.saveCatalog(cat);
372 }
373 
375  : _frameDict(frameDict), _transform(), _pixelOrigin(), _pixelScaleAtOrigin(0 * lsst::geom::radians) {
376  _computeCache();
377 };
378 
379 std::shared_ptr<ast::FrameDict> SkyWcs::_checkFrameDict(ast::FrameDict const& frameDict) const {
380  // Check that each frame is present and has the right type and number of axes
381  std::vector<std::string> const domainNames = {"ACTUAL_PIXELS", "PIXELS", "IWC", "SKY"};
382  for (auto const& domainName : domainNames) {
383  if (frameDict.hasDomain(domainName)) {
384  auto const frame = frameDict.getFrame(domainName, false);
385  if (frame->getNAxes() != 2) {
387  "Frame " + domainName + " has " + std::to_string(frame->getNAxes()) +
388  " axes instead of 2");
389  }
390  auto desiredClassName = domainName == "SKY" ? "SkyFrame" : "Frame";
391  if (frame->getClassName() != desiredClassName) {
393  "Frame " + domainName + " is of type " + frame->getClassName() +
394  " instead of " + desiredClassName);
395  }
396  } else if (domainName != "ACTUAL_PIXELS") {
397  // This is a required frame
399  "No frame with domain " + domainName + " found");
400  }
401  }
402 
403  // The base frame must have domain "PIXELS" or "ACTUAL_PIXELS"
404  auto baseDomain = frameDict.getFrame(ast::FrameSet::BASE, false)->getDomain();
405  if (baseDomain != "ACTUAL_PIXELS" && baseDomain != "PIXELS") {
407  "Base frame has domain " + baseDomain + " instead of PIXELS or ACTUAL_PIXELS");
408  }
409 
410  // The current frame must have domain "SKY"
411  auto currentDomain = frameDict.getFrame(ast::FrameSet::CURRENT, false)->getDomain();
412  if (currentDomain != "SKY") {
414  "Current frame has domain " + currentDomain + " instead of SKY");
415  }
416 
417  return frameDict.copy();
418 }
419 
420 lsst::geom::AffineTransform SkyWcs::_linearizePixelToSky(lsst::geom::Point2D const& pix00,
421  lsst::geom::SpherePoint const& coord,
422  lsst::geom::AngleUnit const& skyUnit) const {
423  // Figure out the (0, 0), (0, 1), and (1, 0) ra/dec coordinates of the corners
424  // of a square drawn in pixel. It'd be better to center the square at sky00,
425  // but that would involve another conversion between sky and pixel coordinates
426  const double side = 1.0; // length of the square's sides in pixels
427  auto const sky00 = coord.getPosition(skyUnit);
428  auto const dsky10 = coord.getTangentPlaneOffset(pixelToSky(pix00 + lsst::geom::Extent2D(side, 0)));
429  auto const dsky01 = coord.getTangentPlaneOffset(pixelToSky(pix00 + lsst::geom::Extent2D(0, side)));
430 
431  Eigen::Matrix2d m;
432  m(0, 0) = dsky10.first.asAngularUnits(skyUnit) / side;
433  m(0, 1) = dsky01.first.asAngularUnits(skyUnit) / side;
434  m(1, 0) = dsky10.second.asAngularUnits(skyUnit) / side;
435  m(1, 1) = dsky01.second.asAngularUnits(skyUnit) / side;
436 
437  Eigen::Vector2d sky00v;
438  sky00v << sky00.getX(), sky00.getY();
439  Eigen::Vector2d pix00v;
440  pix00v << pix00.getX(), pix00.getY();
441  // return lsst::geom::AffineTransform(m, lsst::geom::Extent2D(sky00v - m * pix00v));
442  return lsst::geom::AffineTransform(m, (sky00v - m * pix00v));
443 }
444 
445 lsst::geom::AffineTransform SkyWcs::_linearizeSkyToPixel(lsst::geom::Point2D const& pix00,
446  lsst::geom::SpherePoint const& coord,
447  lsst::geom::AngleUnit const& skyUnit) const {
448  lsst::geom::AffineTransform inverse = _linearizePixelToSky(pix00, coord, skyUnit);
449  return inverse.inverted();
450 }
451 
452 std::shared_ptr<SkyWcs> makeFlippedWcs(SkyWcs const& wcs, bool flipLR, bool flipTB,
453  lsst::geom::Point2D const& center) {
454  double const dx = 1000; // any "reasonable" number of pixels will do
455  std::vector<double> inLL = {center[0] - dx, center[1] - dx};
456  std::vector<double> inUR = {center[0] + dx, center[1] + dx};
457  std::vector<double> outLL(inLL);
458  std::vector<double> outUR(inUR);
459  if (flipLR) {
460  outLL[0] = inUR[0];
461  outUR[0] = inLL[0];
462  }
463  if (flipTB) {
464  outLL[1] = inUR[1];
465  outUR[1] = inLL[1];
466  }
467  auto const flipPix = TransformPoint2ToPoint2(ast::WinMap(inLL, inUR, outLL, outUR));
468  return makeModifiedWcs(flipPix, wcs, true);
469 }
470 
472  bool modifyActualPixels) {
473  auto const pixelMapping = pixelTransform.getMapping();
474  auto oldFrameDict = wcs.getFrameDict();
475  bool const hasActualPixels = oldFrameDict->hasDomain("ACTUAL_PIXELS");
476  auto const pixelFrame = oldFrameDict->getFrame("PIXELS", false);
477  auto const iwcFrame = oldFrameDict->getFrame("IWC", false);
478  auto const skyFrame = oldFrameDict->getFrame("SKY", false);
479  auto const oldPixelToIwc = oldFrameDict->getMapping("PIXELS", "IWC");
480  auto const iwcToSky = oldFrameDict->getMapping("IWC", "SKY");
481 
482  std::shared_ptr<ast::FrameDict> newFrameDict;
483  std::shared_ptr<ast::Mapping> newPixelToIwc;
484  if (hasActualPixels) {
485  auto const actualPixelFrame = oldFrameDict->getFrame("ACTUAL_PIXELS", false);
486  auto const oldActualPixelToPixels = oldFrameDict->getMapping("ACTUAL_PIXELS", "PIXELS");
487  std::shared_ptr<ast::Mapping> newActualPixelsToPixels;
488  if (modifyActualPixels) {
489  newActualPixelsToPixels = pixelMapping->then(*oldActualPixelToPixels).simplified();
490  newPixelToIwc = oldPixelToIwc;
491  } else {
492  newActualPixelsToPixels = oldActualPixelToPixels;
493  newPixelToIwc = pixelMapping->then(*oldPixelToIwc).simplified();
494  }
495  newFrameDict =
496  std::make_shared<ast::FrameDict>(*actualPixelFrame, *newActualPixelsToPixels, *pixelFrame);
497  newFrameDict->addFrame("PIXELS", *newPixelToIwc, *iwcFrame);
498  } else {
499  newPixelToIwc = pixelMapping->then(*oldPixelToIwc).simplified();
500  newFrameDict = std::make_shared<ast::FrameDict>(*pixelFrame, *newPixelToIwc, *iwcFrame);
501  }
502  newFrameDict->addFrame("IWC", *iwcToSky, *skyFrame);
503  return std::make_shared<SkyWcs>(*newFrameDict);
504 }
505 
507  return std::make_shared<SkyWcs>(metadata, strip);
508 }
509 
511  Eigen::Matrix2d const& cdMatrix, std::string const& projection) {
512  auto metadata = makeSimpleWcsMetadata(crpix, crval, cdMatrix, projection);
513  return std::make_shared<SkyWcs>(*metadata);
514 }
515 
517  lsst::geom::Angle const& orientation, bool flipX,
518  lsst::geom::SpherePoint const& boresight, std::string const& projection) {
519  auto frameDict = makeSkyWcsFrameDict(pixelsToFieldAngle, orientation, flipX, boresight, projection);
520  return std::make_shared<SkyWcs>(frameDict);
521 }
522 
524  Eigen::Matrix2d const& cdMatrix, Eigen::MatrixXd const& sipA,
525  Eigen::MatrixXd const& sipB) {
526  auto metadata = makeTanSipMetadata(crpix, crval, cdMatrix, sipA, sipB);
527  return std::make_shared<SkyWcs>(*metadata);
528 }
529 
531  Eigen::Matrix2d const& cdMatrix, Eigen::MatrixXd const& sipA,
532  Eigen::MatrixXd const& sipB, Eigen::MatrixXd const& sipAp,
533  Eigen::MatrixXd const& sipBp) {
534  auto metadata = makeTanSipMetadata(crpix, crval, cdMatrix, sipA, sipB, sipAp, sipBp);
535  return std::make_shared<SkyWcs>(*metadata);
536 }
537 
539  bool simplify) {
540  auto iwcToSky = wcs.getFrameDict()->getMapping("IWC", "SKY");
541  return std::make_shared<TransformPoint2ToSpherePoint>(*iwcToSky, simplify);
542 }
543 
545  auto pixelToIwc = wcs.getFrameDict()->getMapping(ast::FrameSet::BASE, "IWC");
546  return std::make_shared<TransformPoint2ToPoint2>(*pixelToIwc, simplify);
547 }
548 
550  os << wcs.toString();
551  return os;
552 };
553 
554 } // namespace geom
555 } // namespace afw
556 } // namespace lsst
table::Key< table::Array< std::uint8_t > > frameSet
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
static std::shared_ptr< T > dynamicCast(std::shared_ptr< Persistable > const &ptr)
Dynamically cast a shared_ptr.
Definition: Persistable.cc:18
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
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
constexpr double asRadians() const noexcept
Return an Angle&#39;s value in radians.
Definition: Angle.h:166
MatrixMap is a form of Mapping which performs a general linear transformation.
Definition: MatrixMap.h:42
constexpr double asDegrees() const noexcept
Return an Angle&#39;s value in degrees.
Definition: Angle.h:169
lsst::geom::Point2D skyToPixel(lsst::geom::SpherePoint const &sky) const
Compute pixel position(s) from sky position(s)
Definition: SkyWcs.h:349
An object passed to Persistable::write to allow it to persist itself.
An affine coordinate transformation consisting of a linear transformation and an offset.
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< const ast::Mapping > getMapping() const
Get the contained mapping.
Definition: Transform.h:135
bool isFits() const
Return true getFitsMetadata(true) will succeed, false if not.
Definition: SkyWcs.cc:249
table::PointKey< double > crpix
Definition: OldWcs.cc:131
Interface supporting iteration over heterogenous containers.
Definition: Storable.h:58
Point2D getPosition(AngleUnit unit) const noexcept
Return longitude, latitude as a Point2D object.
Definition: SpherePoint.cc:129
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:109
T to_string(T... args)
Channel provides input/output of AST objects.
Definition: Channel.h:60
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< 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:471
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:364
table::PointKey< double > crval
Definition: OldWcs.cc:130
std::shared_ptr< ast::FrameDict > readLsstSkyWcs(daf::base::PropertySet &metadata, bool strip=true)
Read an LSST celestial WCS FrameDict from a FITS header.
ItemVariant const * other
Definition: Schema.cc:56
lsst::geom::SpherePoint getSkyOrigin() const
Get the sky origin, the celestial fiducial point.
Definition: SkyWcs.cc:187
static int constexpr CURRENT
index of current frame
Definition: FrameSet.h:105
bool equals(typehandling::Storable const &other) const noexcept override
Compare this object to another Storable.
Definition: SkyWcs.cc:358
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
STL class.
std::shared_ptr< const ast::FrameDict > getFrameDict() const
Get the contained FrameDict.
Definition: SkyWcs.cc:247
AngleUnit constexpr radians
constant with units of radians
Definition: Angle.h:108
A class representing an angle.
Definition: Angle.h:127
A specialized form of Channel which reads and writes FITS header cards.
Definition: FitsChan.h:202
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:71
Transform< Point2Endpoint, Point2Endpoint > TransformPoint2ToPoint2
Definition: Transform.h:300
STL class.
T sin(T... args)
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
Definition: SkyWcs.cc:362
A WinMap is a linear Mapping which transforms a rectangular window in one coordinate system into a si...
Definition: WinMap.h:45
Eigen::Matrix2d getCdMatrix() const
Get the 2x2 CD matrix at the pixel origin.
Definition: SkyWcs.cc:204
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
AngleUnit constexpr degrees
constant with units of degrees
Definition: Angle.h:109
SkyFrame is a specialised form of Frame which describes celestial longitude/latitude coordinate syste...
Definition: SkyFrame.h:66
A base class for image defects.
bool hasFitsApproximation() const
Does this SkyWcs have an approximate SkyWcs that can be represented as standard FITS WCS...
Definition: SkyWcs.h:364
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:452
Frame is used to represent a coordinate system.
Definition: Frame.h:157
std::shared_ptr< RecordT > src
Definition: Match.cc:48
char delimiter
Definition: Schema.cc:193
std::string toString() const override
Create a string representation of this object.
Definition: SkyWcs.cc:337
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
std::shared_ptr< FrameDict > copy() const
Return a deep copy of this object.
Definition: FrameDict.h:123
static std::shared_ptr< SkyWcs > readStream(std::istream &is)
Deserialize a SkyWcs from an input stream.
Definition: SkyWcs.cc:289
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:260
T str(T... args)
std::shared_ptr< daf::base::PropertyList > getPropertyListFromFitsChan(ast::FitsChan &fitsChan)
Copy values from an AST FitsChan into a PropertyList.
T cos(T... args)
T dynamic_pointer_cast(T... args)
static std::string getShortClassName()
Definition: SkyWcs.cc:279
std::shared_ptr< TransformPoint2ToPoint2 > getPixelToIntermediateWorldCoords(SkyWcs const &wcs, bool simplify=true)
Return a transform from pixel coordinates to intermediate world coordinates.
Definition: SkyWcs.cc:544
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:523
std::shared_ptr< TransformPoint2ToSpherePoint > getIntermediateWorldCoordsToSky(SkyWcs const &wcs, bool simplify=true)
Return a transform from intermediate world coordinates to sky.
Definition: SkyWcs.cc:538
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
constexpr double asArcseconds() const noexcept
Return an Angle&#39;s value in arcseconds.
Definition: Angle.h:178
ShiftMap is a linear Mapping which shifts each axis by a specified constant value.
Definition: ShiftMap.h:40
std::shared_ptr< typehandling::Storable > cloneStorable() const override
Create a new SkyWcs that is a copy of this one.
Definition: SkyWcs.cc:333
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
std::ostream & operator<<(std::ostream &os, Storable const &storable)
Output operator for Storable.
Definition: Storable.h:174
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:198
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
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
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
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
STL class.
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
static std::shared_ptr< SkyWcs > readString(std::string &str)
Deserialize a SkyWcs from a string, using the same format as readStream.
Definition: SkyWcs.cc:317
lsst::geom::Point2D getPixelOrigin() const
Get the pixel origin, in pixels, using the LSST convention.
Definition: SkyWcs.h:215
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
Definition: SkyWcs.cc:366
T pow(T... args)
SkyWcs(SkyWcs const &)=default
Class for storing generic metadata.
Definition: PropertySet.h:67
std::string writeString() const
Serialize this SkyWcs to a string, using the same format as writeStream.
Definition: SkyWcs.cc:327
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition: SkyWcs.cc:506
std::shared_ptr< daf::base::PropertyList > getFitsMetadata(bool precise=false) const
Return the WCS as FITS WCS metadata.
Definition: SkyWcs.cc:218
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
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
String-based source and sink for channels.
Definition: Stream.h:180
int m
Definition: SpanSet.cc:49
lsst::geom::SpherePoint SpherePoint
Definition: misc.h:35
table::PointKey< int > pixel
void writeStream(std::ostream &os) const
Serialize this SkyWcs to an output stream.
Definition: SkyWcs.cc:322
static int constexpr BASE
index of base frame
Definition: FrameSet.h:104
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
table::Key< int > version
Definition: PhotoCalib.cc:374
bool isFlipped() const
Does the WCS follow the convention of North=Up, East=Left?
Definition: SkyWcs.cc:281
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
Extent< double, 2 > Extent2D
Definition: Extent.h:400
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:269
table::Schema schema
Definition: SkyWcs.cc:70
Reports errors from accepting an object of an unexpected or inappropriate type.
Definition: Runtime.h:167
lsst::geom::Angle getPixelScale() const
Get the pixel scale at the pixel origin.
Definition: SkyWcs.h:208
A FrameSet whose frames can be referenced by domain name.
Definition: FrameDict.h:67
std::ostream * os
Definition: Schema.cc:746
STL class.
bool operator==(SkyWcs const &other) const
Equality is based on the string representations being equal.
Definition: SkyWcs.cc:161
AffineTransform const inverted() const
Return the inverse transform.
std::shared_ptr< Object > read()
Read an object from a channel.
Definition: Channel.cc:52
A stream for ast::Channel.
Definition: Stream.h:41
A FrameSet consists of a set of one or more Frames (which describe coordinate systems), connected together by Mappings (which describe how the coordinate systems are inter-related).
Definition: FrameSet.h:99
double constexpr PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:39
bool strip
Definition: fits.cc:901
A class used to convert scalar POD types such as double to Angle.
Definition: Angle.h:70
Transform LSST spatial data, such as lsst::geom::Point2D and lsst::geom::SpherePoint, using an AST mapping.
Definition: Transform.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
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
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
static bool singleClassEquals(T const &lhs, Storable const &rhs)
Test if a Storable is of a particular class and equal to another object.
Definition: Storable.h:151