LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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"
41#include "lsst/pex/exceptions.h"
43
44namespace lsst {
45namespace afw {
46
49
50namespace geom {
51namespace {
52
53int 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
61double const TIGHT_FITS_TOL = 0.0001;
62
63class SkyWcsPersistenceHelper {
64public:
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
81private:
82 SkyWcsPersistenceHelper()
83 : schema(),
84 wcs(schema.addField<table::Array<std::uint8_t>>("wcs", "wcs string representation", "")) {}
85};
86
87class SkyWcsFactory : public table::io::PersistableFactory {
88public:
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
103std::string getSkyWcsPersistenceName() { return "SkyWcs"; }
104
105SkyWcsFactory registration(getSkyWcsPersistenceName());
106
107ast::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
133Eigen::Matrix2d makeCdMatrix(lsst::geom::Angle const& scale, lsst::geom::Angle const& orientation,
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
146std::shared_ptr<TransformPoint2ToPoint2> makeWcsPairTransform(SkyWcs const& src, SkyWcs const& dst) {
147 auto const dstInverse = dst.getTransform()->inverted();
148 return src.getTransform()->then(*dstInverse);
149}
150
151SkyWcs::SkyWcs(daf::base::PropertySet& metadata, bool strip)
152 : SkyWcs(detail::readLsstSkyWcs(metadata, strip)) {}
153
154SkyWcs::SkyWcs(ast::FrameDict const& frameDict) : SkyWcs(_checkFrameDict(frameDict)) {}
155
156bool 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
193Eigen::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
199Eigen::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
264bool 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
295
296bool 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 }
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
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
373bool SkyWcs::equals(typehandling::Storable const& other) const noexcept {
374 return singleClassEquals(*this, other);
375}
376
377std::string SkyWcs::getPersistenceName() const { return getSkyWcsPersistenceName(); }
378
379std::string SkyWcs::getPythonModule() const { return "lsst.afw.geom"; }
380
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
394std::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
435lsst::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
460lsst::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
467std::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
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::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
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
table::Schema schema
Definition python.h:134
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< 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< 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< 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 int constexpr CURRENT
index of current frame
Definition FrameSet.h:105
static int constexpr 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
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
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
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
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:489
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.
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:128
constexpr double asArcseconds() const noexcept
Return an Angle's value in arcseconds.
Definition Angle.h:185
A class used to convert scalar POD types such as double to Angle.
Definition Angle.h:71
Point in an unspecified spherical coordinate system.
Definition SpherePoint.h:57
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:930
ndarray::Array< std::uint8_t, 1, 1 > stringToBytes(std::string const &str)
Encode a std::string as a vector of uint8.
Definition Utils.cc:43
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
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
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
AngleUnit constexpr degrees
constant with units of degrees
Definition Angle.h:110
Extent< double, 2 > Extent2D
Definition Extent.h:400
AngleUnit constexpr radians
constant with units of radians
Definition Angle.h:109
double constexpr PI
The ratio of a circle's circumference to diameter.
Definition Angle.h:40
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